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ABSTRACT 

We have updated and extended our semi-analytic galaxy formation modelling capa- 
bilities and applied them simultaneously to the stored halo/subhalo merger trees of 
the Millennium and Millcnnium-II simulations. These differ by a factor of 125 in mass 
resolution, allowing explicit testing of resolution effects on predicted galaxy proper- 
ties. We have revised the treatments of the transition between the rapid infall and 
cooling flow regimes of gas accretion, of the sizes of bulges and of gaseous and stellar 
disks, of supernova feedback, of the transition between central and satellite status as 
galaxies fall into larger systems, and of gas and star stripping once they become satel- 
lites. Plausible values of efficiency and scaling parameters yield an excellent fit not 
only to the observed abundance of low-redshift galaxies over 5 orders of magnitude 
in stellar mass and 9 magnitudes in luminosity, but also to the observed abundance 
of Milky Way satellites. This suggests that reionisation effects may not be needed to 
solve the "missing satellite" problem except, perhaps, for the faintest objects. The 
same model matches the observed large-scale clustering of galaxies as a function of 
stellar mass and colour. The fit remains excellent down to ~ 30 kpc for massive galax- 
ies. For M* < 6 x 10 10 M Q , however, the model overpredicts clustering at scales below 
~ 1 Mpc, suggesting that the assumed fluctuation amplitude, as = 0.9, is too high. 
The observed difference in clustering between active and passive galaxies is matched 
quite well for all masses. Galaxy distributions within rich clusters agree between the 
simulations and match those observed, but only if galaxies without dark matter sub- 
halos (so-called orphans) are included. Even at MS-II resolution, schemes which assign 
galaxies only to resolved dark matter subhalos cannot match observed clusters. Our 
model predicts a larger passive fraction among low-mass galaxies than is observed, as 
well as an overabundance of ~ 1O 1O M galaxies beyond z ~ 0.6. (The abundance of 
~ 10 11 M Q galaxies is matched out to z ~ 3.) These discrepancies appear to reflect 
deficiencies in the way star-formation rates are modelled. 

Key words: cosmology: theory - cosmology: dark matter mass function - galaxies: 
luminosity function, stellar mass function - galaxies: haloes - cosmology: large-scale 
structure of Universe 



The ACDM model has been successful in interpreting a 
wide variety of observations. These include the cosmic mi- 
crowave background fluctuations at z ~ 1000 (e.g. Dunkley 
et al. 2009), the large-scale clustering of galaxies in the low- 
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1 INTRODUCTION 



redshift universe (e.g. Percival 2010), the cosmic shear field 
measured by weak gravitational lensing (e.g. Fu 2008), the 
high- redshift power spectrum probed by the Lyman a forest 
(e.g. Viel et al. 2009; Paschos et al. 2009), and the abun- 
dance (e.g Vikhlinin et al. 2009) and baryon fractions (e.g 
Allen et al. 2008) of galaxy clusters. Current N-body simula- 
tions can follow the growth of representative samples of dark 
matter halos at high resolution and in their full cosmologi- 
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cal context on scales ranging from those of rich clusters to 
those of dwarf galaxies. The formation of galaxies does not, 
however, trace that of their dark matter halos in a simple 
manner, and the exponentially growing body of high-quality 
galaxy data coming from large surveys cannot be properly 
compared to the ACDM model without a careful treatment 
of baryonic processes. Such detailed comparison is the most 
promising route to clarifying the complex astrophysics un- 
derlying galaxy formation, and it may also uncover prob- 
lems with the ACDM model which are not evident on larger 
scales. 

In the standard scenario of galaxy formation, originally 
proposed by White & Rees (1978), gas cools and condenses 
at the centres of a population of hierarchically merging dark 
matter haloes. These and earlier ideas (e.g. Rees & Ostriker 
1977) were adapted by Blumenthal et al. (1984) to the spe- 
cific initial conditions predicted in a CDM dominated uni- 
verse, showing that the dichotomy between rapid and slow 
cooling regimes provides a natural explanation for the di- 
chotomy between individual galaxies and larger systems like 
groups and clusters. Within such models, and in particu- 
lar within the current standard ACDM model, dark mat- 
ter halos grow through accretion and merging to produce 
a present-day halo mass function which has a very different 
shape from the observed luminosity function of galaxies (e.g. 
Benson et al. 2003). If one nevertheless matches the two, as- 
suming bigger galaxies to live in bigger halos, the ratio of 
halo mass to central galaxy light is found to minimize for 
galaxies similar to the Milky Way and to increase rapidly for 
both more massive and less massive systems. The maximal 
efficiency for converting available baryons into stars is about 
20%, and much lower efficiencies are found for the halos of 
rich clusters or dwarf galaxies (Moster et al. 2010; Guo et al. 
2010). 

A popular explanation for the low efficiency of galaxy 
formation in massive halos is that a supermassive black hole 
releases vast amounts of energy when it accretes gas from 
its surroundings, and that this suppresses cooling onto (and 
hence star formation in) the host galaxy (Silk & Rees 1998; 
Birzan et al. 2004; Croton et al. 2006; Bower et al. 2006). 
For low-mass halos, White & Rees (1978) argued that the 
supernova-driven winds invoked by Larson (1974) to ex- 
plain the metallicities of dwarf galaxies might sufficiently 
reduce the efficiency of galaxy formation. Although this has 
been the preferred theoretical explanation ever since, there is 
still no convincing observational evidence that dwarf galaxy 
winds have the properties needed to reproduce the faint end 
of the galaxy luminosity function in a ACDM universe. Cur- 
rent simulations of dwarf galaxy formation, while predicting 
substantial winds, neverthless suppress star formation much 
less effectively than is required (Sawala et al. 2010). The 
UV and X-ray backgrounds heat the intergalactic medium 
and are also thought to affect galaxy formation in small ha- 
los (Doroshkevich et al. 1967; Couchman & Rees 1986; Efs- 
tathiou 1992; Gnedin 2000; Benson et al. 2002; Hoeft et al. 
2006; Okamoto et al. 2008; Hambrick et al. 2009) 

A related issue is the missing satellite problem. Accord- 
ing to the ACDM model, the halo of the Milky Way ac- 
creted many lower mass halos as it grew, many of which 
should have contained small galaxies. Just as low-mass iso- 
lated halos produce too many dwarf field galaxies unless 
their galaxy formation efficiency is extremely low, so these 



accreted halos overpredict the number of dwarf satellites 
around the Milky Way unless their star formation is simi- 
larly suppressed (Kauffmann et al. 1993). Simulations of the 
growth of such halos revealed correspondingly large numbers 
of surviving dark matter subhalos as soon as their resolution 
was high enough (Klypin et al. 1999; Moore et al. 1999), 
and increasing resolution has predicted ever larger numbers 
of ever smaller objects (Diemand et al. 2007; Springel et al. 
2008; Boylan-Kolchin et al. 2009). The number of known 
satellites of the Milky Way has also increased in recent years, 
through effective use of the Sloan Digital Sky Survey (SDSS) 
to detect extremely low luminosity systems (Willman et al. 
2005; Zucker & et al. 2006; Belokurov & et al. 2007; Ko- 
posov et al. 2008), but the apparent discrepancy between 
the predicted and observed numbers has steadily grown. En- 
vironmental effects undoubtedly play a role in determining 
satellite properties: after a galaxy falls into a larger system, 
its gas may be stripped, leading to a rapid decline in star 
formation, dimming its light and reddening its colour. Tidal 
stripping may remove stars or even destroy the satellite al- 
together, contributing gas to the disk of the central galaxy 
and stars to its stellar halo. Nevertheless, since the subhalos 
survive in the simulations, such disruption cannot explain 
the apparent discrepancy. Many low-mass subhalos must be 
dark if the ACDM model is correct. 

An entirely different resolution of these problems could 
lie in a modification of the ACDM model itself. A number of 
authors have suggested that the suppression of small-scale 
structure expected in a warm dark matter model might re- 
duce the abundance of low-mass halos enough to alleviate 
the tension (e.g. Bode et al. 2001; Zavala et al. 2009; Maccio 
et al. 2010). The strongest constraint on this low-mass cut- 
off currently comes from observations of small-scale struc- 
ture in the high-redshift intergalactic medium, as observed 
through the Ly a forest in quasar spectra. These place an 
upper limit on the cut-off wavelength and the correspond- 
ing halo mass, which implies a lower limit on the mass of 
the dark matter particle (Viel et al. 2008; Boyarsky et al. 
2009a, b). Taken at face value, this recent work appears to 
exclude significant warm dark matter effects on any but the 
very faintest galaxies. 

In recent years, the completion of SDSS has allowed a 
determination of the galaxy stellar mass function down to a 
stellar mass of 1O 7 M0. Above about 1O 8 M0 these mass func- 
tions are robust against incompleteness and cosmic variance 
and have very small uncertainties, other than an overall sys- 
tematic coming from the poorly known stellar initial mass 
function (Baldry et al. 2008; Li & White 2009). The large 
sample size makes it possible to retain small mass function 
errors for subsamples split according to additional galaxy 
properties such as colour and environment (Peng & et al. 
2010). This results in a considerable sharpening of the con- 
straints on galaxy formation within the ACDM model (Guo 
et al. 2010; Sawala et al. 2010), making it timely to reassess 
the viability of current models, in particular their ability to 
reproduce the faint end of the galaxy luminosity (or stellar 
mass) function and the faint satellite abundance around the 
Milky Way. 

In this paper we use the combination of the Millennium 
Simulation (MS, Springel et al. 2005) and the Millennium- 
II Simulation (MS-II, Boylan-Kolchin et al. 2009) to ad- 
dress this issue. The latter has 125 times better mass res- 
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olution and 5 times better force resolution than the MS, 
but follows evolution within a box of 5 times smaller side. 
We update our earlier MS-based galaxy formation models 
(Springel et al. 2005; Croton et al. 2006; De Lucia & Blaizot 
2007, hereafter collectively referred to as DLB07) to include 
a better treatment of a number of physical processes, and 
we apply the improved model to both simulations simulta- 
neously. This allows us to test explicitly how limited res- 
olution affects our results. We demonstrate that together, 
the two simulations enable study of the formation, evolu- 
tion and clustering of galaxies ranging from the faint dwarf 
satellites of the Milky Way to the most massive cD galaxies. 
Uncertain astrophysical processes are strongly constrained 
by the precise, low-redshift abundance and clustering data 
provided by the SDSS. Models consistent with these data 
can be tested against other observational data, notably the 
satellite abundance around the Milky Way, but also, for ex- 
ample, the Tully-Fisher relations of giant and dwarf galaxies 
or the properties of high-redshift galaxy populations. 

Previous generations of semi-analytic galaxy formation 
models have been able to reproduce the properties of ob- 
served galaxy populations in ever increasing detail (White & 
Frenk 1991; Kauffmann et al. 1993; Cole et al. 1994; Kauff- 
mann et al. 1999; Somerville & Primack 1999; Cole et al. 
2000; Springel et al. 2001; Hatton et al. 2003; Kang et al. 
2005; Baugh et al. 2005; Croton et al. 2006; Bower et al. 
2006; De Lucia & Blaizot 2007; Somerville et al. 2008; Font 
et al. 2008; Guo & White 2009; Weinmann et al. 2009). The 
DLB07 model was built for the MS simulation and has been 
extensively compared to the abundance, intrinsic properties 
and clustering of galaxies, both in the local universe and at 
high redshift. These comparisons have generally been lim- 
ited to galaxies with stellar masses of at least 1O 9 M0, cor- 
responding approximately to the resolution limit of the MS. 
When the same model is applied to the MS-II, it signifi- 
cantly overpredicts the observed abundance of galaxies near 
this limit and it substantially overpredicts the abundance at 
lower masses (see Fig. 1). The high- mass cut-off is also at 
slightly larger mass than in the new SDSS data, although 
it was consistent with earlier datasets (Croton et al. 2006). 
Clearly, galaxy formation efficiencies are substantially too 
high at low halo mass in the DLB07 model, and slightly too 
high at high halo mass (see also, for example, Fontanot et al. 
2009) 

In the following section, we revisit the DLB07 model, 
improving the treatment of a number of physical processes 
and retuning the uncertain efficiency parameters to obtain 
a better fit to the new SDSS data on abundance and clus- 
tering. In particular, we change the treatments of supernova 
feedback, of the reincorporation of ejected gas, of the sizes 
of galaxies, of the distinction between satellite and central 
galaxies, and of environmental effects on galaxies. Our pa- 
per is organised as follows. In Sec. 2 we briefly describe the 
two N-body simulations on which we implement our galaxy 
formation model. A detailed description of the semi-analytic 
model itself is presented in Sec. 3. In Sec. 4 we compare both 
the abundance and the clustering of galaxies as a function 
of stellar mass, luminosity and colour to low-redshift data 
from the SDSS. We also compare model predictions to the 
observed abundance of satellite galaxies around the Milky 
Way, to the Tully-Fisher relation of isolated galaxies, and to 
the galaxy number density profiles, stellar mass functions, 
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Figure 1. Stellar mass functions predicted by the galaxy for- 
mation model of DLB07. The green curve is the prediction for 
the MS-II and the red curve is that for the MS. Results for the 
two simulations agree well above 10 9 ' 5 Mq, but resolution effects 
cause an underprediction at lower masses in the MS. Black stars 
show the observed function for SDSS/DR7 with error bars in- 
cluding both counting and cosmic variance uncertainties (Li & 
White 2009; Guo et al. 2010). Blue triangles are results for a low- 
redshift sample (0.0033< z <0.05) from SDSS/DR4 taken from 
Baldry et al. (2008); these are corrected for surface- brightness 
incompleteness, but the error bars do not include cosmic vari- 
ance uncertainties. Clearly the model substantially overpredicts 
the abundance of low-mass galaxies and slightly overpredicts the 
masses of high-mass galaxies. 



and intergalactic light fractions of clusters. A final subsec- 
tion focusses on a few model predictions at high redshift. 
Sec. 5 presents a concluding discussion of our results. 



2 N-BODY SIMULATIONS 

We build virtual catalogues of the galaxy population by 
implementing galaxy formation models on the stored out- 
put of two very large cosmological N-body simulations, 
the Millennium Simulation (MS, Springel et al. 2005) and 
the Millennium-II Simulation (MS-II Boylan-Kolchin et al. 
2009). Both simulations assume a ACDM cosmology with 
parameters based on a combined analysis of the 2dFGRS 
(Colless & et al. 2001) and the first-year WMAP data 
(Spergel et al. 2003). The parameters are fi m — 0.25, 
Q b = 0.045, n A = 0.75, n = 1, a 8 = 0.9 and H = 
73 km s _1 Mpc _1 . These cosmological parameters are not 
consistent with more recent analyses of the CMB data (e.g. 
Komatsu et al. 2010) but the relatively small off-sets are 
not significant for most of the issues addressed in this pa- 
per, with the important exception of the small-scale clus- 
tering analysis of section 4.9.) The parameter which devi- 
ates most from recent estimates is as which is quoted as 
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0.809±0.024 for WMAP7 by Komatsu et al. (2010), disgree- 
ing with the simulation value by almost 4a. As shown by 
Angulo & White (2010) simulations based on the MS cos- 
mology can be scaled to represent neighboring cosmologies 
such as WMAP7 to the precision needed for making galaxy 
catalogues. In future work we will use this to show how pre- 
dictions for galaxy properties are affected by the small pa- 
rameter shifts to WMAP7 and other currently allowed cos- 
mologies. 

Both MS and the MS-II trace 2160 3 particles from red- 
shift 127 to the present day. The MS was carried out in 
a periodic box of side 685 Mpc and the MS-II in a box 
of side 137 Mpc. The corresponding particle masses are 
1.18 x 1O 9 M and 9.45 x 1O 6 M , respectively. The small- 
est halos/subhalos we consider contain 20 bound particles, 
and it will turn out that the MS-II has just sufficient reso- 
lution to study dwarf galaxies as faint as those seen around 
the Milky Way. On the other hand, the large volume of the 
MS makes it possible to study rare objects like rich clus- 
ters and bright quasar hosts. In addition, a comparison of 
the two simulations where both have good statistics allows 
us to study how the limited resolution of the MS affects its 
model galaxy populations. 

The particle data were stored at 64 and 68 times for 
the MS and the MS-II, respectively, with the last 60 being 
identical in the two simulations. At each output time, the 
post-processing pipeline produced a friends-of-friends (FOF) 
catalogue by linking particles with separation less than 0.2 of 
the mean value (Davis et al. 1985). The SUBFIND algorithm 
(Springel et al. 2001) was then applied to each FOF group to 
identify all its bound substructures (subhalos). The merger 
trees which are the basis for our galaxy formation mod- 
elling are then constructed by linking each subhalo found 
in a given output to one and only one descendent at the 
following output (Springel et al. 2005; De Lucia & Blaizot 
2007; Boylan-Kolchin et al. 2009). Note that all our galaxy 
formation models are thus based on the growth and merging 
of the population of subhalos, not on the growth and merging 
of the population of halos. This is an important distinction 
which allows us to build much more realistic models for the 
galaxy population, in particular for its merging rates and 
its clustering, than would otherwise be the case. We refer 
readers to Springel et al. (2005) and Boylan-Kolchin et al. 
(2009) for full descriptions of the two simulations. 

The most massive self-bound subhalo in a FOF group 
is referred to as its main subhalo (sometimes the main halo) 
and usually contains most of its mass. Other subhalos of 
the FOF group are referred to as satellite subhalos. After 
implementation of the galaxy formation model, each FOF 
group hosts a "central galaxy", which sits at the potential 
minimum of the main subhalo. Other associated galaxies 
may sit at the potential minima of smaller subhalos, or may 
no longer correspond to a resolved dark matter substruc- 
ture ("orphans"). These galaxies are collectively referred to 
as satellites, although we note that in this paper we break 
with our previous practice by assuming that the physical 
processes affecting satellite galaxies only begin to differ from 
those affecting central galaxies when a satellite first enters 
the virial radius of the larger system. This is to account for 
the fact that FOF groups quite often link two essentially in- 
dependent dark matter clumps, and the two central galaxies 



are expected to keep evolving quasi-independently while this 
is the case. 

We define the centre of a FOF group to be its potential 
minimum and its virial radius to be the radius of the largest 
sphere with this centre and a mean overdensity exceeding 
200 times the critical value. The total mass within the virial 
radius is defined as the virial mass of the group. Virial radius 
and virial mass are then related by 

1/3 



_/?vir — 



G M vir 
100 H 2 (z) 



(1) 



The virial radius usually lies almost entirely within the 
boundary of the FOF group and, as a result, the virial mass 
is typically somewhat smaller than the FOF mass (and also 
typically somewhat larger than the mass of the main sub- 
halo) . 



3 GALAXY FORMATION MODELS 

Galaxies form at the centres of dark matter halos and gain 
stars by formation from their interstellar medium (ISM) and 
by accretion of satellite galaxies. We assume the ISM to 
form a disk and to be replenished both by diffuse infall from 
the surroundings and by gas from accreted satellite galax- 
ies. Diffuse infall can occur directly from the intergalactic 
medium (in a so-called "cold flow") or through cooling of 
a hot halo surrounding the galaxy. The interaction of these 
processes with each other and with flows driven by super- 
novae and by active galactic nuclei is responsible for the 
overall evolution of each galaxy, which thus cannot be fol- 
lowed realistically without superposing a complex network 
of baryonic astrophysics on the assembly history of its dark 
matter component. Physical understanding of most of these 
baryonic processes is quite incomplete and is based largely 
on simplified numerical simulations and on the phenomenol- 
ogy of observed systems. Descriptions of the processes are 
thus necessarily both approximate and uncertain, and mod- 
els of the kind we build here may offer the best means to 
constrain them empirically using observational data. 

Here we implement simplified galaxy formation recipes 
onto the subhalo merger trees built from the MS and the 
MS-II. Treating baryonic evolution by post-processing cos- 
mological N-body simulations in this way makes it possible 
to explore a wide model and parameter space in a relatively 
short amount of time. In general, our models build on those 
developed in Springel et al. (2005), Croton et al. (2006), 
and De Lucia & Blaizot (2007) hereafter collectively referred 
to as DLB07. The baryonic content of galaxies is split into 
five components, a stellar bulge, a stellar disk, a gas disk, 
a hot gas halo, and an ejecta reservoir. These components 
exchange material through a variety of processes and their 
total mass grows through accretion from the intergalactic 
medium. As noted above the main modifications here con- 
cern the definition of satellite galaxies, a mass-dependent 
model for supernova feedback, the gradual stripping and 
disruption of satellite galaxies, more realistic treatments of 
the growth of gaseous and stellar disks, a model to calcu- 
late bulge end elliptical galaxy sizes, and an updated reion- 
ization model. We determine the free parameters of these 
models using the observed abundance, structure and clus- 
tering of low redshift galaxies as a function of stellar mass, 
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luminosity and colour. In the following we describe our new 
galaxy formation model in detail. For a more general re- 
view of semi-analytic models, we refer the reader to Baugh 
(2006); Benson & Bower (2010b); Benson (2010) 

3.1 Reionization 

It is now well established that the global baryon to dark 
matter mass ratio is 15-20%. In galaxy clusters, the observed 
baryon fraction is close to but somewhat below this value, 
and is mostly in the form of hot gas. In halos like that of the 
Milky Way, on the other hand, at most about 20% of the 
expected baryons are detected and these are primarily in the 
form of stars; the detected baryon fraction is apparently even 
lower in the halos of dwarf galaxies (e.g. Guo et al. 2010). 
One mechanism which may contribute to the low efficiency 
of dwarf galaxy formation is photo-heating of pregalactic 
gas by the UV background. This inhibits gas condensation 
within dark matter halos if the thermal energy exceeds the 
halo potential well depth. This effect was first pointed out 
by Doroshkevich et al. (1967) and was later investigated in 
the context of CDM models by Couchman & Rees (1986) 
and Efstathiou (1992). 

In recent years a number of simulations of this process 
have been carried out. Here we use a fitting function of the 
form originally proposed by Gnedin (2000) to describe how 
the baryon fraction in a halo depends on mass and redshift: 

f b (z,M viI ) = fr i i + (2 q/3 -i) - u 



M F (z) 



(2) 



In this formula, /™ s =17% is the universal baryon fraction as 
given by first-year WMAP estimates (Spergel et al. 2003), 
and a = 2 is a fit to the simulations in Okamoto et al. 
(2008). Mf is the characteristic halo mass of this "filter". 
In halos with M v i r 3> Mf(z) the baryon fraction is set to 
the universal value, while in halos with M v i r -C Mf(z) it 
drops as {M^/Mf) 3 ■ The redshift dependence of Mf is 
determined by the details of how the reionization process 
occurred. Here we use a table of Mf(z) kindly provided by 
Okamoto et al. (2008) from their simulations; Mf varies 
from ~ 6.5 x 10 9 M Q at z = to ~ 10 7 M Q just after reion- 
ization at z ~ 8. These results are consistent with those 
found earlier by Hoeft et al. (2006). In DLB07, a lower 
value of a and a different Mf{z) were adopted (following 
Kravtsov et al. 2004) leading to the substantially higher 
value Mf ~ 3 x 10 10 M Q at z = 0. These differences in sim- 
ulation results appear to reflect differences in resolution and 
in the treatment of radiative transfer. Although we adopt 
the more recent results as "standard" , we will rcdiscuss how 
these issues affect dwarf galaxy formation in Sec. 4.8, show- 
ing that reionization does not appear to be a major factor 
except, possibly, for the faintest Milky Way satellite systems. 



3.2 Cooling 

The pressure of the intergalactic medium has little effect 
on the growth of more massive halos. A fraction ~ f^° s 
of the infalling material is expected to be diffuse gas, and 
thus must shock as it joins the halo. At early times and in 
low-mass halos, post-shock cooling is rapid and the accre- 
tion shock is very close to the central object, which thus 



gains new material at essentially the free-fall rate; at late 
times and in massive halos, post-shock cooling times ex- 
ceed halo sound crossing times, the accretion shock moves 
away from the galaxy, and infalling gas forms a quasi-static 
hot atmosphere from which it condenses onto the central 
galaxy through a cooling flow (Forcada-Miro & White 1997; 
Birnboim &: Dekel 2003). The critical mass separating these 
two regimes is around 1O 12 M0 and is weakly dependent on 
redshift but strongly dependent on f^° s and on the metal- 
licity of the infalling gas (Rees & Ostriker 1977; White & 
Frenk 1991). These rapid infall and quasi-static cooling flow 
regimes have featured in almost all galaxy formation mod- 
els of the last two decades (e.g. Kauffmann et al. 1993; Cole 
et al. 1994; Kauffmann et al. 1999; Somerville & Primack 
1999; Cole et al. 2000; Springel et al. 2001; Hatton et al. 
2003; Kang et al. 2005; Somerville et al. 2008, and also, 
of course DLB07) as well as (by construction) in all direct 
simulations of the galaxy formation process (e.g. Navarro 
& Steinmetz 1997; Steinmetz 1999; Springel & Hernquist 
2003). The simple criterion of White & Frenk (1991) is used 
to separate the two regimes in most semi-analytic models, 
and tests with both one-dimensional (Forcada-Miro & White 
1997) and three-dimensional (Benson et al. 2001; Yoshida 
et al. 2002; Cattaneo et al. 2007) simulations have shown 
it to provide an adequate description. More recent numeri- 
cal work has focussed on the fact that the two regimes can 
effectively coexist near the transition, with cold gas falling 
in narrow streams through a hot gas atmosphere or even a 
galactic wind (see, for example, Fig. 2 of Springel & Hern- 
quist (2003) or Dekel et al. (2009)). A recent reanalysis by 
Benson &: Bower (2010a) emphasised that the details of how 
such "cold flows" are treated has little effect on predicted 
galaxy properties once the necessary strong feedback is in- 
cluded. 

Here we use the simple model of Springel et al. (2001) 
to estimate the gas cooling rate in the hot halo regime. We 
assume that infalling gas is shock-heated to the virial tem- 
perature of the host halo at an accretion shock, and that its 
distribution interior to this shock is a quasistatic isothermal 
sphere with density falling as the inverse square of radius. 
The cooling time at each radius can then be calculated as 



tcooi(r) = 



3fimnkT v i 



(3) 



2phot(r)A(T hot ,Z hot )' 

where /jiiih is the mean particle mass, k is the Boltz- 
mann constant, phot{r) is the hot gas density at radius r, 
A(Thot, Zhct) is the temperature- and metallicity-dependent 
cooling function (Sutherland & Dopita 1993), and Zhot is the 
metallicity of the hot halo gas. Tkot = 35.9(VW/km s _1 ) 2 K 
is the assumed virial temperature of the host halo. For main 
subhalos, the gas temperature is updated according to the 
current circular velocity at the virial radius at each snapshot, 
while for satellite subhalos, we assume the gas temperature 
to be constant at the value it had when the subhalo was 
accreted onto its main halo. 

The cooling radius is then estimated through 



T'cool 



tdyn,hWlhotA(Tli 0t , ^hot) 



6nfj,m}ikT v i T R v 



(4) 



The definition of the halo dynamical time id yn ,h involves an 
arbitrary constant. Here we adopt the convention id yn ,h = 
-Rvir/Kir = 0.1H(zy 1 as in De Lucia et al. (2004). Readers 
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are refereed to Croton et al. (2006) and Somerville et al. 
(2008) for the discussion of other possible choices of idyn.h 
when defining r CO oi- When r coo i < Rvir, we assume that we 
are indeed in the cooling flow regime, and that the cooling 
rate onto the central galaxy is 1 



Mcooi = nihot 



^"cool 



Rvir ^dyn,h 



(5) 



When r coo i > Rvir, on the other hand, we assume that 
we are in the rapid infall regime and gas accretes onto the 
central object in free-fall, thus on the halo dynamical time: 



m-hot 



dyn,h 



(6) 



Note that condensation is smoother in time in this model 
(which is essentially identical to that of De Lucia et al. 
(2004)) than in the model of DLB07, who assumed the 
hot gas to fall onto the central object instantaneously as 
soon as it satisfies r coo i > -Rvir- In a situation of steady ac- 
cretion onto a low-mass halo, the DLB07 model results in 
non-steady behaviour; after a cooling "event" empties the 
halo, its hot gas atmosphere is replenished by infall until it 
again reaches the rapid cooling threshold, triggering another 
cooling event. Although the time-average condensation rate 
is equal to the infall rate onto the halo, condensation oc- 
curs in bursts which induce (possibly) unrealistic star for- 
mation bursts in the central object. The model of Equ. (6) 
eliminates this behaviour. For a low-mass halo experiencing 
steady infall, condensation onto the central object is now 
also steady, and the hot gas atmosphere has constant mass 
equal to the gas infall rate times the halo dynamical time. 
The coefficient of unity in Equ. (6) ensures continuity of the 
condensation rate as a halo transitions between the rapid 
infall and hot halo regimes. With these changes, condensa- 
tion rates onto galaxies fluctuate strongly only in response 
to variations in the accretion onto their halos, not as a con- 
sequence of discontinuities in the way we treat the various 
regimes. 

Another substantive difference in the way we treat cool- 
ing with respect to the model of DLB07 is that we now allow 
satellite galaxies to have their own hot gas halos which can 
be removed dynamically by tidal and ram-pressure effects. 
This hot gas can continue to cool onto the (satellite) galaxy, 
adding to its interstellar medium and providing additional 
fuel for star formation. We discuss this in more detail in 
sect. 3.6 below. 



3.3 Disk Sizes 

Disk sizes are not only interesting in their own right, but are 
also important because they determine the surface density 



1 Note that the coefficient on the rhs of this equation differs from 
that in the corresponding equation (equ. 6) of Croton et al. (2006). 
By checking the original code, we have verified that a factor of 
0.5 was erroneously omitted when programming equ. 6 in this 
paper, and that this error then propagated through all the later 
DLB07 papers, with the result that the equations which actually 
correspond to the models of Croton et al. (2006) and DLB07 are 
those given here. For consistency in comparing to the earlier work, 
we have kept these assumptions in our new model. 



of gas in disks, which in turn determines the star forma- 
tion rate. DLB07 followed the simple model of Mo et al. 
(f998) which assumes that the specific angular momentum 
of a galaxy disk is the same as that of the dark halo in which 
it is embedded. This results in the characteristic size of a 
disk scaling as the product of the virial radius and the spin 
parameter of its host halo. Mo et al. (f 998) intended this as 
a simple model for a population of isolated disk galaxies at a 
single time, and several difficulties arise when it is applied to 
individual objects as they grow in time. For example, halo 
spin parameters often change discontinuously by quite large 
amounts as new material is accreted, but it is not plausible 
that this should result in instantaneous changes in size of 
the central disk. Here we introduce a new and more realistic 
disk model which distinguishes between gas and stellar disks 
and allows each of them to grow continuously in mass and 
angular momentum in a physically plausible way. 

We assume that the change in the total vector angu- 
lar momentum of the gas disk during a timestep can be 
expressed as 



AJ E1 



5J, 



gas,coolinf 



+ 5J K . 



<5Jgas,SF, 



(7) 



where <5J gas , cooling, <5J g as,acc and <5Jsf are respectively the 
angular momentum changes due to addition of gas by cool- 
ing, to accretion from minor mergers, and to gas removal 
through star formation. 

When new gas condenses onto the central galaxy, we as- 
sume it to carry specific angular momentum which matches 
the current value for the dark matter halo Jdm/Mdm- The 
angular momentum change due to this gas can thus be ex- 
pressed as 



SJg&s. cooling — -^cool , j St, 

Mom 



(8) 



where M coo i is the condensation rate from Equ. (5) or 
Equ. (6), and 5t is the timestep. When a minor merger hap- 
pens (which we define as the smaller galaxy having a bary- 
onic mass less than a third that of the larger) we assume 
the cold gas from the smaller object to be added to the disk 
of the larger (see Sec. 3.7), carrying specific angular mo- 
mentum equal to the current value for the dark matter halo 
of the larger object. The corresponding angular momentum 
change in the gas disk is thus 



gas,acc — -/^/sa-t,gas , * 
IVh 



(9) 



where M sa t,gas is the cold gas mass in the satellite disk. 

When some gas from the cold disk is converted into 
stars, we assume it to have the average specific angular mo- 
mentum of the gas disk, J sas /M gas , so the change in angular 
momentum of the gas and stellar disks is 



<5Xas,SF = —M* /? as St = — <5J„,SF 
M Eas 



(10) 



where M» is the star formation rate. 

When the cold gas in disks is reheated by SN feedback, 
we assume that the outflowing material also carries away 
its "fair share" of the angular momentum. As a result, the 
specific angular momentum of the gas disk is not changed 
by the SN feedback process. 

For the stellar disk, we assume the total change in (vec- 
tor) angular momentum over the timestep to be oV^sf- Thus 
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we are assuming that the angular momentum of the disk is 
changed only by star formation. In particular, bulge forma- 
tion through disk instability (see Sect. 3.8) is assumed to 
remove negligible angular momentum from the disk. 

We assume both the gas disk and the stellar disk to be 
thin, to be in centrifugal equilibrium and to have exponential 
density profiles 



E(i? g 
and 



E ga sO exp( — -Rgas/-Rgas,d), 



S(i?») = E, exp(-i?»/il J ,,d) 



(11) 



(12) 



where i? gas ,d and i?*,d are the scale-lengths of the gas and 
stellar disks, and E gas o and E*o are the corresponding cen- 
tral surface densities. Assuming a flat circular velocity curve, 
as would hold for a galaxy with negligible self-gravity in an 
isothermal dark matter halo, the scale-lengths can be calcu- 
lated as 

J gas / -A^gas 



R 



gas,d 



2V m 



and 



R*.d — 



J./M,, d 

2 Vmax 



(13) 



(14) 



where M ga s and M,,d are the total masses of the two disks. 
A significant issue here is the dark matter response to the 
gravity of the baryons as they condense within the halo. The 
simplest model for this effect assumes adiabatic contraction 
within a spherical halo (e.g. Barnes 1984; Blumenthal et al. 
1986) and has often been adopted in galaxy formation mod- 
els (e.g. Mo et al. 1998; Cole et al. 2000). However, recent 
simulations suggest that this simple scheme overestimates 
the effect (e.g. Gnedin et al. 2004; Abadi et al. 2010). Most 
recently, Tissera et al. (2010) found disk maximum rota- 
tion velocities very similar to the maximum circular veloci- 
ties found for dark matter only haloes formed from identi- 
cal ACDM initial conditions. Here, for simplicity, we adopt 
Vmax, the maximum circular velocity of the surrounding dark 
matter halo, as the typical rotation velocity for both stellar 
and HI disks. Note that we keep the rotation velocities of 
satellite galaxies fixed after infall. This is because the in- 
ner regions of the dark matter subhalo, which determine the 
rotation velocity of the disks, change rather little until the 
subhalo is about to be destroyed (e.g. Hayashi et al. 2003; 
Kazantzidis et al. 2004). 

Fig. 2 shows a few results for this simple model imple- 
mented on the MS-II to demonstrate that it gives results 
in fair agreement with observation. The first panel com- 
pares model predictions for the distribution of stellar half- 
mass radius for disk galaxies as a function of their stellar 
mass to observational results from the Sloan Digital Sky 
Survey (SDSS, Shen et al. 2003). This SDSS study defined 
"late-type" galaxies according to the concentration param- 
eter c = Rgo/Rso, where R$o and Rqo are the radii which 
enclose 50% and 90% of the projected stellar light. Galaxies 
with c < 2.86 were taken to be late-type, primarily spiral 
galaxies. To calculate R50 and -R90 for our model galaxies, 
we assume the above exponential model for the disk and 
a Jaffe profile for the bulge (the modelling of bulge size 
will be described in Sec. 3.8). In practice, c < 2.86 corre- 
sponds roughly to galaxies with M^d/M^tot > 0.80 in our 
model. The solid curves in the figure show the median and 



±lcr range of the model distribution as a function of stellar 
mass, while the symbols show the SDSS data. The ampli- 
tude, slope and scatter of the observations are all fairly well 
reproduced, although the predicted slope is somewhat shal- 
lower than observed. 

The second panel of Fig. 2 shows the distribution of the 
ratio of the sizes of the gas and stellar disks. It gives the 
median, the upper and lower quartiles and the upper and 
lower 10% points as a function of stellar mass for the same 
model galaxies plotted in the first panel of the figure. Gas 
disks are typically larger than stellar disks by about a factor 
of 1.6 but the scatter in the ratio is large. This agrees quite 
well with the observational situation. For the WHISP sample 
of Noordermeer et al. (2005) the average ratio of disk sizes 
is 1.72 and values within their sample of 49 galaxies range 
from 0.6 to 4.1. 

The third panel of Fig. 2 shows the distribution of the 
misalignment angle 8 = arccos (j gas ■ J*/\ Jgas 1 1 J* |) between 
the two disks for several ranges of stellar mass, again for 
the same galaxies. The distribution of misalignment angles 
is quite broad and seems to depend very little on stellar 
mass. Warps are quite often seen in the outer parts of spiral 
galaxies, particularly when the outer HI distribution is com- 
pared to the inner stellar disk. The structure and evolution 
of the two components is quite strongly coupled (e.g. Bin- 
ney & Tremaine (2008) section 6.6), but our simple model 
nevertheless gives some indication of the extent to which 
misalignments might reflect changes with time in the orien- 
tation of accreted angular momentum. 



3.4 Star Formation 

In this paper we assume stars to form from cold gas in the 
disk according to a simplified version of the empirical rela- 
tion which Kennicutt (1998) found to give a good descrip- 
tion of galaxy-scale star formation in the bulk of low-redshift 
star-forming galaxies. Stars form efficiently only in regions 
where the surface mass density exceeds a critical value which 
is plausibly related to the Toomre (1964) threshold for local 
instability of a rotationally supported disk. Toomre's crite- 
rion is a function of local velocity dispersion, of the surface 
densities of stars and gas, and of the local rotation curve. We 
adopt a simple model which assumes a flat rotation curve 
and a gas velocity dispersion which is everywhere 6 km/s, 
leading to the critical surface density suggested in Kauff- 
mann (1996) and Croton et al. (2006), 



E crit (i?) = 12 x 



( Vm 



\ 200km/s J \ lOkpc 



( R 



M pc" 



(15) 



We integrate this out to three exponential scale radii i? gas ,d 
and then divide by a factor of 2 to obtain a critical gas mass 
which is required for any stars to form 



M C1 



11.5 x 10' 



gas.d \ 



^OOkm/s J ylOkpc J 



(16) 



The final reduction by a factor of 2 is introduced to agree 
with the assumptions of Croton et al. (2006) who took the 
cold gas surface density to be constant with radius in disks 
at threshold. 

The amount of cold disk gas that is converted into long- 
lived stars per unit time is assumed to be 
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9.5 



1 1.5 



10.0 10.5 11.0 
log(M.(M sun )) 

Figure 2. The top panel gives the distribution of the radius con- 
taining half the stellar mass as a function of stellar mass for local 
late-type galaxies. These are defined as having SDSS concentra- 
tion parameter c < 2.86 (see details in the text). The solid curve 
is the median half-mass radius predicted by our model applied to 
the MS-II, while dashed curves indicate the rras scatter in log R 
at each stellar mass. Symbols are the observed median and scat- 
ter from the SDSS study by Shcn ct al. (2003). The central panel 
gives the 10, 25, 50, 75 and 90% points of the distribution of 
the ratio of sizes of the gaseous and stellar disks in our model, 
also as a function of total stellar mass, while the bottom panel 
shows the same percentile points of the distribution of the relative 
inclination of the two disks. 



a(M gas - M cri t)/tdyn 



(17) 



where tdyn = S-Rgas.d/Vmax is the characteristic timescale at 
the edge of the star-forming disk, and a is an adjustable 
efficiency parameter. We will adopt a = 0.02, which results 
in a few percent of the gas in a disk being converted into 
stars each rotation period. The star formation rates implied 
by this model are, in the mean, quite similar to those in 
DLB07, but our revised treatments of cooling and of disk 
size lead to considerably smoother evolution than before, 



with less "bursty" star formation histories in the bulk of the 
galaxy population. 



3.5 Supernova Feedback 

During their short lives, massive stars emit large amounts 
of radiation through optical and UV emission, and large 
amounts of mechanical energy through their winds. As they 
die, comparable amounts of radiative and mechanical energy 
are liberated by the final supernova (SN) explosion. This can 
dramatically reshape the surrounding interstellar medium, 
ionising and heating it, and in many cases driving galactic- 
scale winds, Such effects are generically referred to as SN 
feedback. As Larson (1974) showed, they can have a major 
impact on the evolution of low-mass galaxies, determining, 
for example, their metallicities. White & Rees (1978) argued 
that such SN feedback may induce the strong dependence 
of galaxy formation efficiency on halo mass required to ex- 
plain why most stars live in galaxies with stellar mass close 
to the upper limit imposed by cooling constraints, and why 
the overall conversion of baryons into stars is relatively in- 
efficient. These ideas have subsequently been explored by 
many authors, particularly in the context of understanding 
the shape of the galaxy luminosity function at low lumi- 
nosities (e.g. Benson et al. 2003). Here we assume that SN 
feedback injects gas from the cold disk into the hot halo and, 
in addition, can transfer halo gas to the ejecta reservoir. 

We estimate the amount of cold disk gas that is reheated 
by SN feedback and injected into the hot halo component 



<5M rchoat = e disk x SM t . 



(18) 



where 5M* is the mass of newly formed long-lived stars. 
DLB07 took e d isk to be a constant, based on some observa- 
tional indications that mass outflow rates are typically a few 
times the star formation rate in actively star-forming galax- 
ies. We find that this scaling does not suppress star forma- 
tion in low-mass galaxies enough to reproduce the shallow 
slope of the observed stellar mass function, so we have ex- 
tended it to allow higher ejection efficiencies in dwarf galax- 
ies, taking 



Edisk 



e X 



0.5 



+ ( 



V m 



^70km/s 



-Pi 



(19) 



where e and /3i are free parameters describing the ratio of 
reheated mass to new stellar mass in massive galaxies, and 
the scaling of this ratio with Knax in dwarfs. The circular 
velocity dependence is motivated by the fact that less energy 
is needed to heat a solar mass of gas to the halo virial tem- 
perature and to eject it from the disk in lower mass galaxies. 
A naive argument leads to the expectation j3\ ~ 2, but a va- 
riety of factors could lead to a different dependence, so we 
adjust both fi\ and e when fitting to observations, in par- 
ticular to the observed stellar mass function. Below we will 
take e = 6.5 and j3\ = 3.5 in our standard model. 

We parametrise the total amount of energy effectively 
injected by massive stars into disk and halo gas as: 



1 2 
<5-Esn = thaio x -<5M»V SN . 



(20) 



where 0.5V SN is the mean kinetic energy of supernova ejecta 
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per unit mass of stars formed, and, following Croton et al. 
(2006), we take Vsn = 630km/s, based on standard assump- 
tions for the stellar initial mass function and the energetics 
of SN explosions. In this case also, DLB07 assumed the effi- 
ciency e h aio to be a constant. However, since dwarf galaxies 
have stronger winds, lower metallicities and less dust than 
galaxies like our own, it is plausible that radiative losses 
during the thermalisation of ejecta energy are substantially 
smaller in dwarfs than in giants. We have therefore allowed 
for this possibility explicitly by setting 



/ Kir 
inc = —7 I ■= 



ehalo = V X 



0.5 + 



/ Vm 



^ 70km/s 



-Pi 



(21) 



where r\ is a free parameter which encodes possible variations 
about our IMF and SN assumptions, possible energy input 
from the winds and UV radiation of massive stars, and the 
radiative losses during ejecta thermalisation, while /?2 de- 
scribes the dependence of this last factor on Vmax- Again, 
we adjust these parameters when fitting to observations of 
the stellar mass function and gas-to-star ratios of galaxies. 
Our standard model below adopts r) = 0.32 and /?2 = 3.5. 
We expect that ehalo < 1 and, unlike DLB07 or Bower et al. 
(2006), we enforce this constraint in all our models. 

Given this energy input into disk and halo gas, the total 
amount of material that can ejected from a halo/subhalo can 
be estimated as 



5M e jec = 



5-Esn — 5<5M re heatK 2 ir 
2 vir 



(22) 



If this equation gives <5M e j ec < 0, we assume that the mass 
of reheated gas saturates at <5M rchcat = 5Esn/(^V? t ) and 
that no gas is ejected from the halo/subhalo. 

The reheating and ejection efficiencies, edisk and thaio, 
decline with increasing halo circular velocity, saturating at 
0.5e and 0.5r/, respectively. This dependence is controlled 
by the values of j3\ and P2 which are chosen to fit the abun- 
dance of low-mass galaxies. j3i primarily affects the low-mass 
slope of the stellar mass function, while P2 affects its ampli- 
tude. Our default model has a very strong V ma x-dependence, 
/3i = /?2 = 3.5, but because of saturation effects the results 
are not very sensitive to this choice. For example, the adop- 
tion of /3i = /?2 = 1-5 predicts a stellar mass function only 
slightly steeper than our default model and overpredicts the 
abundance of galaxies of stellar mass 1O 8 M by 0.1 dex com- 
pared to the default model. This dependence of SN feedback 
on Vmax also affects the metallicities of low-mass galaxies 
(see details in Sec. 4.3). Compared to DLB07, our model 
gives stronger feedback at low circular velocities. This is the 
primary reason that it produces fewer dwarf galaxies and 
that these have lower metallicities than in the earlier model. 

The gas mass M c j oc thrown out of a system by these 
effects is stored in an ejecta "reservoir" associated with the 
halo/subhalo, whence it may later be reincorporated into the 
hot halo gas and so again become available for cooling onto 
the central galaxy. In low-mass halos, hot winds are likely 
to leave at a substantially higher velocity relative to the 
escape velocity and so their gas is likely to be more difficult 
to reaccrete. To allow for this possibility, we introduce a 
dependence of the reaccretion rate on halo/subhalo virial 
velocity, 



\ 220km/s I Wdynji 



(23) 



where 7 is a free parameter which we take to be 0.3. With 
these choices, ejected gas is returned to the hot halo com- 
ponent in a few dynamical times for galaxies like the Milky 
Way, but takes substantially longer in dwarf systems. 

The association of hot gaseous halos and ejecta reser- 
voirs with satellite subhalos is a substantial change in our 
model with respect to DLB07. As detailed in the next sub- 
section we explicitly model the stripping of these compo- 
nents by tidal and ram-pressure effects. 



3.6 Satellite Galaxies in Groups and Clusters 

In the following, we classify galaxies into three types accord- 
ing to their relationship to the dark matter distribution. 
Type galaxies are the central galaxies of main subhalos 
and so can be considered as the principal galaxies of their 
FOF groups. Type 1 (satellite) galaxies lie at the centre 
of non-dominant subhalos, while type 2 (satellite) galaxies 
are those which no longer have an associated dark matter 
subhalo which is resolved by the simulation. The latter are 
often referred to as "orphan galaxies" . All galaxies are born 
as type 0. They usually become type 1 when they fall into 
a group or cluster, and they may later become type 2. Type 
2's may later merge into the central galaxy of their halo. 

FOF halos often link together two (or more) essen- 
tially disjoint dark matter structures, joining them with low- 
density "bridges" of particles. In such a situation, one would 
expect the central galaxies of the various objects to evolve 
independently until the smaller ones actually fall into the 
main body of the system. To represent this we have changed 
the treatment of type 1 galaxies from that in DLB07. When 
a type galaxy first becomes type 1 (i.e. its FOF halo is first 
linked to a more massive FOF halo) we continue to treat it 
as a type galaxy (i.e. in the same manner as a galaxy 
at the centre of a main subhalo) until it falls within R v i r 
of the centre of its new FOF halo. At this point we switch 
on tidal and ram-pressure stripping processes which can re- 
move gas from the galaxy or even disrupt it completely. In 
our model such processes only occur within _R v i r so that if 
a satellite passes outside R v i T again it is once more treated 
in the same way as a type galaxy. 2 This change primarily 
affects galaxies between i? v i r and the boundary of the FOF 
group. It leads to a reduction in the number of "true" satel- 
lite galaxies (e.g. galaxies whose evolution is influenced by 
being a non-central object within a larger system). We dis- 
cuss the number of galaxies affected by this change, as well 
as the overall number of satellites and of orphans in our MS 
and MS-II models, in Appendix A. 



2 Since only the main subhalo of a FOF halo has an associated 
R v ; r , this quantity is not available for "independent" type 1 galax- 
ies outside the _R v i r of their new FOF halo. For such objects we 
use the values of i? v i r and M v ; r recorded at the last time they 
were type O's when values of these quantities are required in our 
cooling recipe. 
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3.6.1 Gas Stripping 

In most semi-analytic models, hot gas associated with a halo 
is assumed to be stripped immediately after accretion onto 
a larger system, leading to a rapid decline in star formation 
and a reddening in colour (e.g. Wang et al. 2007; Weinmann 
et al. 2006; Baldry et al. 2006; Font et al. 2008). In the 
real Universe (Sun et al. 2007; Jeltema et al. 2008) and in 
hydrodynamic simulations, however, the hot atmosphere of 
massive satellite galaxies may survive for some considerable 
time after accretion. McCarthy et al. (2008) found that for 
satellite galaxies with typical structural and orbital param- 
eters, up to 30% of the initial hot halo gas can remain in 
place for up to 10 Gyr. Weinmann et al. (2009) and Font 
et al. (2008) constructed MS-based models for incremental, 
rather than instantaneous removal of material through tidal 
stripping and ram-pressure stripping. In the following we 
describe how we include both mechanisms in our own mod- 
els, which are similar to but different in detail from those of 
Weinmann et al. (2009) and Font et al. (2008). 

We assume the hot gas in a subhalo to have a distri- 
bution that exactly parallels that of the dark matter. Since 
tidal acceleration acts identically on hot gas and dark mat- 
ter at each location, we take the hot gas mass to be reduced 
by tidal stripping in direct proportion to the subhalo's dark 
matter mass. The latter is, of course, followed explicitly in a 
dynamically consistent way by the original simulation. Thus 
we assume 

Afhot(-Rtidal) Mdm 



Mr 



(24) 



where MDM.infaii and Mhot.infaii were the DM mass of the 
subhalo and the mass of its associated hot gas when its cen- 
tral galaxy was last a type 0, and Mdm and M hot are the 
current masses of these two components. Recall that we as- 
sume p oc r~ 2 for the hot gas distribution, thus M^ ot {r) oc r. 
The tidal radius beyond which the hot gas is stripped can 
be thus expressed as 



Mom 



Mi 



DM, infall 



-Rdm, 



(25) 



where -RDM.infaii was the virial radius of the subhalo just 
before it became a satellite. 

In addition to tidal forces, the hot gas around satellite 
galaxies experiences ram-pressure forces due to satellite's 
motion through the intracluster medium (ICM). At a certain 
distance, R T . P ., from the centre of the satellite, self-gravity 
is approximately balanced by this ram pressure: 

psat (J2r.p.)Kat — Pp ar (fl)Vorbit, (26) 

where p Ba .t{Rr.p.) is the hot gas density of the satellite at 
radius R r .p., V sa ,t is the virial velocity of the subhalo at in- 
fall (which we assume to be constant as the subhalo orbits 
around the main halo), p pa . T (R) is the hot gas density of the 
parent dark matter halo at distance R from the centre of its 
potential well, and Vorbit is the orbital velocity of the satel- 
lite, which we estimate as the virial circular velocity of the 
main halo. The ram-pressure dominates over gravity beyond 
R T .p. and hot gas at these radii is stripped. 

We compare the two radii -Rtidai and i? r . P . and define 
the minimum of the two as the stripping radius 



/ 1 



mm(7?tidai, Rv. P . 



(27) 



Beyond i? s tri P , we assume all the hot gas to be removed with- 
out modifying the gas profile within i? s tri P • Thus the cooling 
rate onto the centre is not affected, to lowest order, by this 
stripping. We assume gas in the "ejecta reservoir" of a sub- 
halo to be stripped in proportion to the hot gas. It is unclear 
where this reservoir should be located and whether or not 
it will be affected by ram-pressure effects (e.g. whether it 
is primarily diffuse or in relatively dense clouds). Thus we 
adopt the simple approach of stripping it in proportion to 
the hot gas. Stripped material from each of these compo- 
nents is added to the corresponding component associated 
with the central (type 0) galaxy of the main subhalo, and 
so can never fall back onto its original subhalo. 

In addition to stripping, at least two other processes 
affect the hot gas halos of satellites. One is cooling. The 
hot gas in satellite galaxies can cool onto the central cold 
star-forming disk. We assume that the temperature of the 
hot gas atmosphere is not changed by stripping and cool- 
ing processes, remaining pegged to its value at infall. The 
cooling rate is calculated just as in Sec. 3.2, which ensures 
continuity in cooling rates as central galaxies turn into satel- 
lite galaxies. As cooling depletes the hot atmosphere, we as- 
sume its density to drop everywhere, but its profile shape 
and bounding radius to remain the same. 

SN feedback also modifies the hot atmospheres and 
ejecta reservoirs of satellite galaxies. As in central galaxies, 
star formation in satellites releases large amounts of energy, 
reheating both cold ISM gas and the hot gas atmosphere. 
Font et al. (2008) presented a model in which both the hot 
and the reheated gas of satellites is stripped primarily in 
the initial infall event. They found that the satellite galaxy 
properties are very sensitive to the way the secondary re- 
heated gas (which is only reheated after the galaxy has be- 
come a satellite) is stripped. If it is stripped as in the initial 
infall event, satellite galaxies lose all gas and become red 
very rapidly. In order to retain gas and keep satellites blue 
for longer, they adopted a stripping efficiency for this sec- 
ondary reheated gas which is only 10% of that at infall, 
and they do not strip any of the hot gas after the initial 
stripping event. In our model, we have adopted continuous 
stripping, in which hot and ejected components are stripped 
equally at each timestep as long as the galaxy is a satellite. 
We assume the reheated gas to extend to a radius equal to 
the virial radius of the subhalo at infall. Taking into account 
the stripping mechanisms discussed above, only reheated gas 
within i?stri P (i.e. a fraction -R s tri P /-RDM,mf of the total re- 
heated gas) remains in the subhalo; the rest is added to the 
hot atmosphere of the main halo. If SN is strong enough 
to predict that material should be ejected from the subhalo 
altogether, then we use the formulae given above for cen- 
tral galaxies (Equations (20) and (22)), and distribute the 
ejected material between the ejecta reservoirs of the satellite 
and central galaxies in the same proportions as the reheated 
gas. In general, the stripping of gas in our model is more effi- 
cient than in Font et al. (2008), but considerably less efficient 
than in DLB07. 

Our current model differs from that of DLB07 both be- 
cause galaxies effectively become satellites later (when they 
cross R v i r rather than when they become part of a larger 
FOF group) and because satellites retain their hot gas com- 
ponents and ejecta reservoirs until these are removed by 
stripping (rather than losing them as soon as they become 
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satellites). Satellites thus retain more fuel for star formation 
and can be expected to stay blue longer. Note that ram- 
pressure stripping does not affect the cold gas component of 
galaxies in our models. This is unrealistic for galaxies in the 
inner regions of rich clusters and results in passive SO galax- 
ies retaining significant gas and dust which should probably 
have been removed (see Fig. 11). 

We illustrate the effect of our modification of stripping 
recipes in Fig. 3. We select 1000 galaxy clusters from the MS 
with Mvir > 2 x 10 14 Mq. For each, we calculate the fraction 
of actively star-forming galaxies as a function of projected 
distance from the centre in units of R v i r , and we average 
over all clusters. "Actively star-forming" here means having 
a specific star formation rate (SSFR, the ratio of star for- 
mation rate to stellar mass) above 10 _11 yr _1 . We consider 
galaxies with velocity relative to cluster centre (peculiar + 
HqX line-of-sight distance difference) less than 3 x V v ir, 
dividing them into four stellar mass bins as indicated by 
the log M* ranges given in the bottom right corner of each 
panel. To emphasize the environmental effects which con- 
cern us here, each active galaxy fraction is normalized to its 
"field" value, estimated at 20i? v ir- Symbols with error bars 
are observational data from Weinmann et al. (2009) based 
on the SDSS cluster sample of von der Linden et al. (2007). 
Predictions from our model are in red, those from the model 
of DLB07 in black. Clearly, within R v i T , the changes we have 
introduced do slow the decline of star formation in satellite 
galaxies, although the differences are not large. This is be- 
cause satellites start to be stripped later in our new recipes, 
and even thereafter they retain some hot gas which can fuel 
star formation and keep them blue, whereas in DLB07, hot 
gas is stripped immediately once galaxies become attached 
to a larger FOF group and star formation ceases once their 
cold ISM gas is used up. Note, however, that the fraction 
of active galaxies in the field differs between our model and 
DLB07, with our model predicting somewhat lower active 
fractions in general, thereby worsening the overall agreement 
with observation. This is a result of the enhanced feedback 
we have introduced in order to match the observed stellar 
mass function (see Sec 3.9). 



3.6.2 Disruption 

The stellar component in subhalos can also be stripped in 
the presence of very strong tidal forces. Usually, the galaxy is 
harder to disrupt than its dark halo because it is more com- 
pact and denser. We thus assume that the stellar component 
of a satellite galaxy is affected by tidal forces only after its 
subhalo has been entirely disrupted, i.e. it has become a 
type 2 galaxy. The position of such a galaxy is identified 
with that of the most bound particle of its subhalo at the 
last time the subhalo could be identified. To estimate when 
stripping of stars is important we assume the satellite orbits 
in a singular isothermal potential, 



4>{R) = Kir hi R. 



(28) 



Assuming conservation of energy and angular momentum 
along the orbit, its pericentric distance can be estimated 
from: 
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Figure 3. The reduction in the fraction of actively star-forming 
galaxies (M*/M* > 10 _11 yr~ 1 ) as a function of projected dis- 
tance from cluster centre in units of ii v i r . The four panels refer to 
different ranges of log M*/Mq as indicated by the labels. Predic- 
tions from the preferred model of this paper applied to the MS are 
shown in red, those from the model of DLB07 in black. Symbols 
with error bars are SDSS data for a large sample of nearby clusters 
taken from Weinmann et al. (2009). For each curve the fraction of 
actively star-forming galaxies is normalised by its "field" value, 
taken to be the value at 20 R v i T . This emphasises the effect of 
cluster environment on star formation activity. 



where R is the current distance of the satellite from halo 
centre, and V and Vt are the velocity of the satellite galaxy 
with respect to halo centre and its tangential part, respec- 
tively. 

We compare the main halo density at pericentre with 
the average baryon mass (cold gas mass + stellar mass) den- 
sity of satellite within its half mass radius. If 
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(30) 



(29) 



we assume the satellite galaxy is disrupted entirely. Its stars 
are then assigned to a population of intracluster stars (ICS) 
and its cold gas and the associated metals are added to the 
hot gas atmosphere of the halo central galaxy. Note that this 
calculation does not fully account for dynamical friction ef- 
fects on the satellite orbit, which are underestimated by the 
simulation once the remaining mass of a subhalo drops below 
the stellar mass of its associated galaxy. Note also that we do 
not model continuous disruption. Rather, once Equ. (30) is 
satisfied, satellite galaxies are disrupted completely. When 
a central type galaxy merges in to a larger system and 
becomes a type 1 satellite, it carries its ICS with it until it 
becomes a type 2 galaxy. At this point, its current central 
galaxy acquires its ICS. 
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3.7 Mergers 

Mergers can occur between a central galaxy and a satellite 
galaxy, and between two satellite galaxies. In the MS, the 
minimum resolved subhalo has a mass of 2.3 x 1O 1O M . The 
stellar mass of the galaxy within a given subhalo is thus 
smaller than the subhalo mass, except in the case of very 
massive satellites. In the MS-II, however, the minimum sub- 
halo mass is 1.9 x 10 s M , and the stellar mass of a galaxy 
often becomes larger than the mass of its host subhalo well 
before we lose the track of the subhalo. In this situation 
the simulation no longer correctly follows the expected de- 
cay of the satellite orbit through dynamical friction. In this 
paper we therefore modify the DLB07 treatment of merg- 
ers, which estimated a dynamical friction time until merging 
only once the satellite's subhalo is fully disrupted. Here we 
estimate this decay time as soon as the mass of a subhalo 
drops below that of the galaxy it contains, and we immedi- 
ately set the countdown clock for merging. The position and 
velocity of the satellite galaxy are thereafter traced by the 
most bound particle of the subhalo at the time the merger 
clock was switched on, modified by a time-dependent orbit- 
shrinking factor which models the orbital decay caused by 
the dynamical friction (see below). As in DLB07, we adopt 
the dynamical friction formula of Binney & Tremaine (1987) 
to estimate the merging time for a satellite galaxy, 



^friction — Ctfri 



1C Gm sat In A ' 



(31) 



where afH C = 2.34. DLB07 found this coefficient to be 
needed to reproduce observed luminosity functions at the 
luminous end. It was also found to be appropriate in the 
N-body studies of Boylan-Kolchin et al. (2008) and Jiang 
et al. (2008). Unlike DLB07, we here take m aat to be the 
sum of the baryonic mass of the satellite galaxy and the 
dark matter mass of its subhalo. The dark matter mass here 
refers to the mass of the subhalo just before the merger clock 
is set. r sat is the distance between the central and satellite 
galaxies at the time when we start the merger clock, and 
In A = ln(l + M v i r /m sat ) is the Coulomb logarithm. After a 
time teiction the satellite galaxy is assumed to merge with 
the central galaxy of the main halo. If a main halo is ac- 
creted onto a larger system and becomes a subhalo, any of 
its satellites for which the merger clock is already set are 
assumed to keep orbiting within this subhalo and to merge 
into its central galaxy when the time runs out. In this way, 
our model allows satellites to merge into the central galax- 
ies of both main and subdominant subhalos (although the 
latter is quite rare). 

We have also attempted to model approximately the dy- 
namical friction driven orbital decay of type 2 galaxies which 
leads to their eventual merging with the central galaxy, 
even though the low-mass subhalo or the simulation par- 
ticle with which the galaxy is associated clearly experiences 
no such decay. A simple model in which an "isothermal" 
satellite spirals to the centre of a larger "isothermal" host 
on a near circular orbit predicts that the radius of the orbit 
should decay linearly in time. To mimic this, we multiply 
the positional offset of the tracer particle from the central 
galaxy with which its galaxy is destined to merge by a factor 
(1 — At/tfriction) to define the position of the galaxy, where 
At is the time since the merger clock was initialised. The 



velocity of the galaxy is kept equal to that of the tracer 
particle. 

Our modelling differentiates between major and minor 
mergers. Major mergers are those between galaxies with 
baryonic masses differing by less than a factor of three. More 
extreme mass ratios are treated as minor mergers. During 
a major merger, the disks of the progenitors are destroyed 
completely, leading to the formation of a spheroidal rem- 
nant. In a minor merger, the disk of the larger progenitor 
survives and accretes the cold gas and stellar components of 
the small galaxy. In both cases, the merger triggers a star- 
burst which we represent using the "collisional starburst" 
model of Somerville et al. (2001). During the merger, a frac- 
tion, e b urst, of the cold gas of the merging galaxies is con- 
verted into stars, where 



eburst = 0.56 



M min 
M mai 



(32) 



and M m i nor and M ma j or are the total baryon masses of the 
minor and major progenitors, respectively. The coefficient 
and index here are consistent with those given by Cox et al. 
(2008) and Somerville et al. (2008). The stars formed during 
major mergers contribute to the elliptical remenants, while 
those formed during minor mergers are added to the disks. 
Feedback and chemical enrichment from the starburst are 
modeled in the same way as for quiescent star formation, 
and the strong SN feedback produced by a major merger 
can expel almost all the remaining cold gas from the sys- 
tem, suppressing further star formation until a new gas disk 
grows. 

To summarize, our treatment of mergers differs from 
that of DLB07 only in that we switch on the merger clock 
as soon as the dark matter mass of a subhalos drops below 
the baryonic mass of its central galaxy, that we take into 
account the baryonic mass of the galaxy when calculating 
the dynamical friction time, and that we include an approx- 
imate representation of the shrinkage of orbits caused by 
dynamical friction. Many aspects of these recipes are quite 
crude but they nevertheless represent reasonably well the 
results of recent simulations of both gas-poor and gas-rich 
galaxy mergers (e.g. Naab & Burkert 2003; Cox et al. 2008). 



3.8 Bulge Formation 

Three modes of bulge growth are included in our model: 
major mergers, minor mergers and disk buckling. 

After a major merger, all stars from the progenitors 
and all the newly formed stars are assumed to end up in 
a spheroidal component. After a minor merger, the disk of 
the larger progenitor remains intact but its bulge acquires all 
the pre-existing stars from the minor progenitor, while the 
newly formed stars are added to the disk. In both cases, the 
spheroidal component grows in mass and changes in size. We 
use energy conservation and the virial theorem to estimate 
the change in size: 
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(33) 



where C is a structure parameter relating the binding energy 
of a galaxy to its mass and radius, and center is a parame- 
ter quantifying the effective interaction energy deposited in 
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Figure 4. The distribution of galaxies across morphological type 
as a function of stellar mass. Red lines show the fraction of galax- 



ies with 



A/to 



> 0.7, which we consider to represent elliptical 



galaxies. Blue lines indicate galaxies with 0.03 < 



(normal spirals) and green indicate 



M 



< 0.03, essentially 



A/total 

pure-disk or extreme late-type galaxies. Model results for the MS 
are shown with dashed lines and for the MS-II with solid lines. 
The symbols give observational results for real galaxies from Con- 
selice (2006). 
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Figure 5. Half-mass radius as a function of stellar mass for early- 
type galaxies, which we define as galaxies with SDSS concentra- 
tion parameter c > 2.86. The solid curve gives the median half- 
mass radius predicted by our model at each stellar mass, while 
dashed curves show the rms scatter in logi?. Symbols with error 
bars indicate the median and rms scatter of observational esti- 
mates taken from the SDSS study of Shen et al. (2003). 



the stellar components. C = 0.49 for an exponential disk 
whereas C = 0.45 for a bulge with an r 1 ' 4 profile; to sim- 
plify, we adopt C = 0.5. We also set center = 0.5, so that 
Qinter/C = 1. This is roughly consistent with the numerical 
simulation results of Boylan-Kolchin et al. (2005) which give 
1.3 < aintcr/C < 1.7 for the most probable orbits of dissi- 



pationless major mergers of elliptical galaxies. We prefer a 
slightly smaller value of center because it gives bulge sizes 
in better agreement with observation (see below). A simi- 
lar formula was used by Cole et al. (2000), but rather than 
following subhalo mass directly, as in our approach, they 
used analytic formulae to estimate the dark matter mass 
of satellites just prior to merging. This dark matter mass 
was included in their analogue of Equ. 33. In our model, the 
dark matter mass of satellites is almost always very small 
(or zero) at the time they merge. Further, we assume the 
final merger to be from a tightly bound orbit. Thus we ne- 
glect the effects of dark matter and use the stellar masses of 
the two objects in Equ. 33. 

The term on the left-hand side of Equ. (33) is the bind- 
ing energy of the final bulge: M ncWi buig C is its stellar mass 
and i?now, bulge is its half-stellar-mass radius. The first and 
second terms on the right-hand side are the self-binding en- 
ergies of the two individual progenitors, while the third term 
is the binding energy invested in their orbit at merger. For 
major mergers, M\ and M2 are the sum of the mass of stars 
and of the cold gas converted into stars for the two progen- 
itors, and Ri and R2 are the corresponding half mass radii. 
For minor mergers, Afi and Ri are the mass and half-mass 
radius of the bulge of the major progenitor, and M2 and R2 
are the stellar mass and the half-stellar-mass radius of the 
minor progenitor. 

Secular evolution is thought to be another important 
channel for the formation of galaxy bulges, in particular in 
systems where the self-gravity of the disk is dominant. Here 
we adopt the same simple, schematic criterion as DLB07 to 
delineate disk instability: 



V m 



< 



GMt, 



(34) 



where M* ; d and R* t d are the stellar mass and exponential 
scale length of the stellar disk, and Vmax, as usual, is the 
maximum circular velocity of the DM (sub)halo hosting the 
disk. In the original presentation of this criterion by Efs- 
tathiou et al. (1982) the factor of 3 was missing and V ma x 
represented the maximum circular velocity of the combined 
disk-halo system. The smaller coefficient used here reflects 
the facts that this latter Vmax is expected to be significantly 
larger than the maximum circular velocity of the unper- 
turbed dark halo for realistic systems near the instability 
boundary, and that more recent simulations have shown ex- 
ponential disks in NFW halos to be somewhat more stable 
than would be inferred from the early experiments of Efs- 
tathiou et al. (1982) (see, for example, Sellwood & Moore 
1999; Sellwood & Evans 2001). 

When the criterion of Equ. (34) is met, we transfer 
mass, 5M*, from the disk to the bulge to keep the disk 
marginally stable. Recall that we assume an exponential 
density profile for the stellar disk. Here we further assume 
that the mass is transferred from the inner part of the disk 
and that the bulge formed in this way occupies the corre- 
sponding region (i.e. the bulge half-mass radius equals to 
the radius of this region): 

&M« = 27rE» i?,,d[i?,,d - (Rb + R*, d )exp(-R h /R*, d )}, (35) 

where Rb is the half-mass radius of the newly formed bulge, 
and covers the region from which the stellar mass is trans- 
ferred into the bulge. We assume that negligible angular 
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momentum is transferred to the bulge from the disk with 
these stars so that the angular momentum of the disk is 
unchanged. Since we also assume an unchanged rotation ve- 
locity and an exponential profile, the disk exponential scale 
length increases and its central surface density decreases 
when stars are transferred to the bulge. 

If there is already a bulge present when a disk goes 
unstable, we assume the instability to produce a new bulge 
with half mass radius Rh given by Equ. (35), which "merges" 
into the existing bulge in the same way as in galaxy mergers, 
simply replacing Mi and Ri with the mass and half-mass 
radius of the existing bulge, and replacing M2 and R2 with 
8M* and R^. The only difference is that we set Qmtcr = 2 in 
this case, since the interaction energy between the "old" and 
"new" bulges is stronger than in the case of galaxy mergers 
since the two are concentric. 

To illustrate how well these recipes work, Fig. 4 com- 
pares observational data to model predictions for the dis- 
tribution of galaxies across morphological type as a func- 
tion of stellar mass. Red curves are for galaxies with 
AfBuigo/Af tota i ^ 0.7 ("elliptical galaxies"), blue for galax- 
ies with 0.7 > Meuigo/Mtotai ^ 0.03 ("normal spirals") and 
green for galaxies with Meuige/Mtotai < 0.03 ("pure disks"). 
Solid and dashed curves are results based on the MS-II and 
the MS, respectively. The two simulations produce conver- 
gent results above log M* = 10, but they differ progressively 
at lower stellar masses because the resolution of the MS 
is no longer good enough to follow accurately the detailed 
formation histories of the galaxies. The symbols in Fig. 4 
are observational results from Conselice (2006). These agree 
well with the models, provided the MS-II results are taken 
at low stellar masses. To study the relative roles in of disk 
instability and mergers in building bulges, we calculated a 
model without the disk instability mode. This showed that 
in our default model, disk instability is a major contribu- 
tor to bulge formation in intermediate mass galaxies like 
the Milky Way. At both higher and lower masses, mergers 
are the dominant mechanism; in particular, massive ellipti- 
cal galaxies are built by mergers. This is consistent with the 
results found by De Lucia et al. (2006) using our previous 
galaxy formation model. 

To illustrate how well our simple recipe reproduces the 
sizes of the spheroidal components of galaxies, Fig. 5 plots 
half-mass radius against stellar mass for early-type galaxies 
defined to be those with concentration parameter c > 2.86 
(see 3.3 for how we estimate c; in practice, this limit cor- 
responds approximately to M Bu i E e/-Mt ota i > 0.20 ). A solid 
curve gives our model prediction for the median half-mass 
radius at each stellar mass, while dashed curves indicate 
the predicted scatter. The symbols are SDSS results for the 
median and scatter from Shen et al. (2003). Overall, agree- 
ment is fair, at best. At lower masses, our default model 
predicts a larger median value than is observed, perhaps 
reflecting gas dissipation during gas-rich mergers (e.g. Hop- 
kins et al. 2009). Indeed, recent work suggests that includ- 
ing gas dissipation may exlain both the steep slope of the 
size vs. stellar mass relation (Dekel & Cox 2006), and its 
small scatter (Covington et al. 2008). It is also noticeable 
that the scatter in size is larger in our simple model than 
in the SDSS data, particularly at low mass. These deficien- 
cies actually become worse if we restrict the sample to more 
strongly bulge-dominated galaxies, since a significant part 



of the trend in this figure is due to the size-stellar mass re- 
lation for disks highlighted in Fig. 2. The small observed 
scatter has recently been confirmed for a large sample of 
visually classified SDSS galaxies by Nair et al. (2010) who 
emphasise that a tight relation may be difficult to under- 
stand if spheroids are built stochastically through mergers. 
Our model confirms that this may be a problem and that a 
more detailed theoretical treatment is warranted. 

3.9 Black Hole Growth and AGN feedback 

There is growing evidence that galactic nuclear activity is 
closely related to galaxy formation. Here we follow Croton 
et al. (2006), separating black hole growth into two modes: 
"quasar" mode and "radio" mode. 

The quasar mode applies to black hole growth during 
gas-rich mergers. During a merger, the central black hole of 
the major progenitor grows both by absorbing the central 
black hole of the minor progenitor, and by accreting cold 
gas. The total growth in mass is calculated as 

SM.u = (tj) ( TTT5 ^L_) ,,36) 

where MBH.min is the mass of the black hole in the minor 
progenitor, M co id is the total cold gas in the two progeni- 
tors, and M m i n and M ma j are the total baryon masses of the 
minor and major progenitors, respectively. Here / is a free 
parameter, which, following Croton et al. (2006) we set to 
0.03 in order to reproduce the observed local Mbh — Mb„i ge 
relation. Both major mergers and gas-rich minor mergers 
contribute significantly to the growth of the black hole. We 
do not explicitly model feedback due to this mode, which al- 
ways coincides with a starburst in the merging galaxies. Any 
feedback from accretion onto the black holes can be thought 
of as being part of the substantial energy input which we 
assume this starburst to produce. As noted above, this is 
often sufficient to eject all the gas from the merger remnant. 

Radio mode growth occurs through hot gas accretion 
onto central black holes. The growth rate in this mode is 
calculated, following Croton et al. (2006), as 

where, for a main subhalo, the hot gas fraction, /hot, is the 
ratio of hot gas mass, Mhot, to subhalo dark matter mass 
Mom, while for a type 1 galaxy in a satellite subhalo, / ho t is 
the ratio of hot gas mass to dark matter mass within i? s t r ip, 
p i?3trlp — M D m infaii- The parameter k sets the efficiency of 

^"-DM,infaIl ' 

hot gas accretion. Again following Croton et al. (2006), we 
assume that this hot mode accretion deposits energy in rela- 
tivistic jets with 10% efficiency and that this energy is then 
deposited as heat in the hot gas atmosphere, as is observed 
directly through radio bubbles in galaxy clusters (e.g. Mc- 
Namara & Nulsen 2007; Birzan et al. 2004). Specifically, we 
assume an energy input rate: 

Sradio = 0.1M B HC 2 , (38) 

where c is the speed of light. The effective (net) mass cooling 
rate is thus 

A/c 00 i jCf f = M coo i — 2i5 rad i /V v 2 ir . (39) 
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Figure 6. The relation between black hole mass and bulge stellar 
mass at z = 0. Red contours give predictions from our model 
applied to the MS-II. The distributions in black hole mass are 
normalised to unity at each stellar mass and the contours indicate 
their 5, 10, 25, 75, 90 and 95 percentage points. The green curve 
represents the median values. Blue crosses are observational data 
taken from Haring & Rix (2004). 



Note this specific form of "radio mode feedback" is only one 
possible description of AGN feedback, and other forms can 
give quite similar results (e.g. Bower et al. 2006). Moreover, 
there is no direct observational evidence that radio mode 
feedback can offset cooling in halos with mass as low as 
~ 1O 13 M0, where it plays an important role in our models. 

In our preferred model the accretion efficiency k is set 
to be 1.5 x 10 -5 in order to match the high-mass end of 
the stellar mass function. This is twice the value adopted in 
DLB07 (k = 7.5 x 10~ 6 ). There are three reasons for this 
change, in addition to the fact that the new SDSS stellar 
mass functions cut off at high mass more strongly than the 
data used in DLB07. The first is that DLB07 assumed the 
hot gas mass of a halo to be f^° a M v i r minus the baryonic 
masses of all the galaxies associated with the FOF group, 
even those which lie outside R v i r . Here we substract only the 
baryonic masses of the galaxies that lie inside 7? v ir, resulting 
in higher estimates of A/hot and so larger cooling rates which 
the radio mode feedback must offset. The second reason is 
that we have introduced a "disruption" mechanism which 
destroys some type 2 galaxies which previously survived in 
galaxy clusters. The ISM of these disrupted galaxies adds 
additional metal-rich material to the hot gas atmosphere, 
again enhancing its predicted cooling rate relative to the pre- 
vious model. The final reason is that the enhanced feedback 
at low mass, which we have introduced in order to match 
the observed z = stellar mass function, results in more 
gas remaining available to cool at later times. Note that our 
model assumes the hot gas in all systems to be distributed 
with p oc r~ 2 at the virial temperature T V1T . In reality, feed- 
back both from star formation and from an AGN may well 
change the profile of the surrounding hot gas, making it less 
centrally concentrated and less able to cool. This would re- 
sult in less need for feedback at later times. (See Bower et al. 
(2008) for a simple model based on this idea.) As may be 



seen in Fig. 6, the increased feedback efficiency in our new 
model does not significantly affect its fit to the observed re- 
lation between the black hole mass and bulge stellar mass. 
This is because black hole growth is in any case dominated 
by the quasar mode. 

Radio mode feedback works in essentially the same way 
in our model as in Croton et al. (2006) and DLB07. It is 
more effective at low redshift and in massive objects, both 
because the black hole is more massive, and because the hot 
gas fraction is higher there. The effect has a very weak, if any, 
dependence on large-scale environment (Croton & Farrar 
2008) . Note that our model differs from DLB07 in that radio 
mode can also operate in satellite galaxies at the centres of 
their own subhalos. In DLB07, such satellite subhalos no 
longer retained any hot gas so that radio mode activity was 
completely quenched there. 

3.10 Metal Enrichment 

Our treatment of metal enrichment follows that of De Lucia 
et al. (2004) quite closely. Here we briefly summarise the 
various processes we include. As stars evolve, both heavy 
elements and a fraction of the initial stellar mass are re- 
turned instantaneously to the cold gas component of the 
ISM. The new material is assumed to be fully mixed with 
the pre-existing cold gas. A more realistic treatment should 
take into account the time delay between star formation 
and the return of both mass and metals to the interstel- 
lar medium. While the return of mass and metals from SN 
type II is indeed effectively instantaneous for the purposes 
of galaxy evolution, the same is not true for SN type la. In 
addition, mass loss and metal enrichment from intermedi- 
ate mass stars takes place over Gyr timescales and is also 
important for a detailed understanding of metallicity pat- 
terns in galaxies. We intend to implement these processes in 
future work. In our current model, metals are carried into 
hot gas atmospheres and ejecta reservoirs when SN feedback 
reheats cold disk gas and ejects it. Metals from both these 
components can then be stripped from satellite galaxies and 
added to the corresponding components of the host system. 
Reincorporation and cooling can then take the metals into 
another (or the same) galaxy again. A more detailed descrip- 
tion of metal enrichment and the exchange between different 
components can be found in De Lucia et al. (2004). 

3.11 Stellar Population Synthesis and Dust 
Extinction 

To compare model prediction with observations, we need 
to calculate the photometric properties of our model galax- 
ies. Here we again follow DLB07, using stellar population 
synthesis models from Bruzual & Chariot (2003). We adopt 
a Chabrier initial function which has fewer low-mass stars 
than a Salpeter IMF and is a better fit to observational 
data both in our own Galaxy and in those nearby early-type 
galaxies for which detailed dynamical data are available (e.g. 
Cappellari et al. 2006). A detailed description can be found 
in De Lucia et al. (2004) . We also follow DLB07 and adopt 
a slab dust model to account for the extinction of the star 
light. At higher redshift, we extend this model as in Guo & 
White (2009) . Extinction is modeled as a function of gas col- 
umn density, metallicity and redshift. The main difference 
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from DLB07, is that a redshift dependence is introduced so 
that for galaxies of given gas metallicity, the gas-to-dust ra- 
tio is higher at high redshift than in the local universe. This 
is motivated by observational data on high-redshift galax- 
ies (e.g. Steidel et al. 2004; Quadri et al. 2008). Kitzbichler 
& White (2007) found that such a redshift dependence was 
needed for the DLB07 galaxy formation model to reproduce 
faint galaxy counts and redshift distributions, while Guo & 
White (2009) needed it to reproduce the abundance and 
clustering pf colour-selected galaxy populations at redshifts 
of 2 and 3 (see these papers for details) . 



4 SYSTEMATIC PROPERTIES OF THE 
GALAXY POPULATION 

In the last section we set out our new galaxy formation 
model and clarified the areas where it significantly alters or 
extends the earlier model of DLB07. Several of these exten- 
sions involve processes which were not previously included, 
notably the separate evolution of the sizes and orientations 
of gaseous and stellar disks, the size evolution of spheroids, 
tidal and ram-pressure stripping of satellite galaxies, and 
the disruption of galaxies to produce intracluster light. To 
illustrate the effects of these new ingredients, we have al- 
ready presented a number of results from a simultaneous 
application of the new model to the MS and MS-II. In the 
current section we present a wide range of further results, 
primarily for the low-redshift universe where recent data 
now constrain the galaxy population over a range exceeding 
six orders of magnitude in stellar mass. By combining the 
MS and MS-II we are able to test our model against ob- 
servation over this full range - the first time this has been 
possible using a direct simulation technique. By combining 
the two simulations we are also able to check explicitly how 
our results are effected by their limited mass resolution (as 
already done, for example, in Fig. 4). 

We begin with a comparison of our model with the ob- 
served stellar mass function of galaxies, because we use this 
as the primary constraint on the various parameters in our 
star-formation and feedback models. We summarize these 
parameters and the values we assign to them in our pre- 
ferred model in Table 1. Note that other model parame- 
ters (for example, those in our treatments of cooling, of disk 
and spheroid sizes, and of stripping, merging and disrup- 
tion) also affect the stellar mass function, but we have set 
these to agree with other simulation or observational data, 
whereas the parameters in Table 1 were chosen primarily to 
fit the stellar mass function, and secondarily to ensure that 
gas-to-star ratios are in reasonable accord with observation. 
Because of the coupling between different parts of the model, 
an iterative method has to be used to find acceptable param- 
eter sets. Those of our preferred model are almost certainly 
not unique, but they all lie within the physically plausible 
range discussed, for example, by Croton et al. (2006). In- 
deed, where the meanings correspond, our parameters are 
close to those presented in that paper and DLB07, except 
in a few cases which we highlight individually. 








log 10 (M.[M o ]) 

Figure 7. The abundance of galaxies as a function of their stellar 
mass. In the upper panel, green and red curves give the predic- 
tions of our preferred model when applied to the MS-II and the 
MS, respectively. The error bars on the MS-II curve show a "cos- 
mic variance" uncertainty estimated from the rms scatter in the 
mass functions among 125 disjoint subvolumes of the MS, each 
with volume equal to that of the MS-II. Stars with error bars are 
the observational result, including cosmic variance uncertainties, 
for SDSS/DR7 as given by Li & White (2009) after a correction 
to total stellar masses following Guo et al. (2010). Blue triangles 
with error bars are the SDSS/DR4 results of Baldry et al. (2008). 
These are corrected for surface brightness incompleteness, but the 
error bars do not include cosmic variance uncertainties which are 
quite large for these low-mass objects. In the lower panel, black 
(Li & White 2009) and blue (Baldry et al. 2008) symbols show the 
abundance ratio of the SDSS data to our model prediction based 
on the MS-II. The red curve is the ratio of our MS and MS-II 
predictions, while the purple curve is the ratio of the DLB07 pre- 
diction to our MS-II prediction. A dashed green line indicates a 
ratio of unity. 



4.1 Stellar Mass and Luminosity Functions 

In Fig. 7 we compare the predictions of our preferred model 
to the observed abundance of galaxies as a function of stellar 
mass. The solid green curve is the prediction based on the 
MS-II, while the solid red curve is based on the MS. The two 
converge well above a stellar mass of about 3 x 1O 9 M0, but 
at lower masses the MS underpredicts abundances because 
it does not resolve halos less massive than 2.3 x M 10 M Q , a.s 
compared to 1.9 x 10 8 in the MS-II. At the highest masses, 
the two simulations also diverge, but this is mainly due to 
cosmic variance and the relatively small volume of the MS- 
II. We estimate this uncertainty by dividing the MS into 
125 sub-cubes, each of the same volume as the MS-II. The 
rms scatter among their individual stellar mass functions is 
given by the error bars overplotted on the green MS-II curve. 
Boylan-Kolchin et al. (2009) show that such differences be- 
come more prominent at high redshift. Black stars are the 
observed stellar mass function estimated from SDSS/DR7 
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Table 1. Summary of those parameters of our preferred model which were adjusted primarily to fit the observed 2 = stellar mass 
function. 



Parameter Description Preferred value 

a Star formation efficiency (Sec. 3. 4) 0.02 

e Amplitude of SN reheating efficiency (Sec. 3.5) 6.5 

01 Slope of SN reheating efficiency (Sec. 3.5) 3.5 
r) Amplitude of SN ejection efficiency (Sec. 3.5) 0.32 

02 Slope of SN ejection efficiency (Sec. 3.5) 3.5 
7 Ejecta reincorporation efficiency (Sec. 3.5) 0.3 

K Hot gas accretion efficiency onto black holes (Sec. 3.9) 1.5 X 10~ 5 



by Li & White (2009), except that the masses are converted 
to total stellar masses as described in Appendix A of Guo 
et al. (2010). Note that these observed stellar masses are 
also based on the same Chabrier IMF used in our models, 
so that IMF uncertainties should not affect the comparison 
of the two. The error bars here include cosmic variance un- 
certainties and are very small, reflecting the large volume of 
the survey. Blue triangles are estimates based on SDSS/DR4 
from Baldry et al. (2008). These include a correction for sur- 
face brightness incompleteness, which becomes significant at 
these low masses, but their error bars do not include cosmic 
variance which is quite large because of the small volume ef- 
fectively surveyed for such faint galaxies. To make the com- 
parison clearer, the lower panel of Fig. 7 shows the SDSS 
data (the symbols), our MS prediction (red curve) and the 
MS predictions of DLB07 (purple curve) all ratioed to our 
predictions based on the MS-II. 

It is clear from Fig. 7 that adopting the MS-II stellar 
mass function below about 3 x 1O 9 M0 and the MS function 
at higher masses results in a very good match to the obser- 
vational results for our preferred parameters. The fit extends 
over the full range of the observations from about 1O 12 M0 all 
the way down to about 2 x 1O 7 M0. The slope at the low-mass 
end is around -1.46 in the model, significantly steeper than 
the value of -1.155 quoted by Li & White (2009). The Baldry 
et al. (2008) results suggest that this may reflect the onset 
of incompleteness effects at the lowest masses considered by 
Li & White (2009). The high resolution of the MS-II allows 
us to predict galaxy abundances to substantially lower stel- 
lar masses. Here we show our predictions down to 1O 6 M 
although there are currently no reliable observations with 
which to compare them. At this mass, the predicted num- 
ber density is around 0.3 Mpc _3 (log M*) _1 . Galaxies even 
less massive than this can be observed in the Local Group 
and we show below that our model does, in fact, agree quite 
well with the abundance of Milky Way satellites as a func- 
tion of luminosity (see Sec. 4.8). 

At high stellar masses, where growth is limited by AGN 
feedback as in Croton et al. (2006), our model overpredicts 
the abundance found by Li & White (2009). This likely re- 
flects the observational difficulty in estimating stellar masses 
for the most luminous cD galaxies in clusters. As a result 
of the problems with dealing with extended envelopes and 
crowded fields, SDSS photometry gives luminosities for such 
galaxies which are significantly lower (by up to one magni- 
tude) than found in other investigations (e.g. von der Linden 
et al. 2007). As the lower panel of Fig. 7 shows, our model 
agrees with these SDSS data significantly better than the 
older model of DLB07. 



In Fig. 8 we show predictions of this same preferred 
model for galaxy luminosity functions in the SDSS g, r, i 
and z bands, comparing them with observational data from 
a low-redshift SDSS sample taken from Blanton et al. (2005). 
In all these plots we have used results from MS+MS-II at ab- 
solute magnitudes brighter than -20, where results from the 
two simulations converge, and results from the MS-II alone 
at fainter magnitudes. Given that Fig. 7 shows our model to 
overpredict slightly the abundance of low-mass dwarf galax- 
ies, it is somewhat surprising that it turns out to underpre- 
dict their abundance as a function of luminosity in all four 
bands. Several effects may contribute to this discrepancy. 
One is that, as we will see later, the fraction of non-star- 
forming dwarf galaxies appears to be significantly larger in 
the model than in the SDSS data, so we are probably assign- 
ing stellar mass-to-light ratios which are too large to many 
dwarfs. A second is that Blanton et al. (2005) corrected 
their luminosity functions for incompleteness at low surface- 
brightness, and their corrections may be larger than those 
applied by Baldry et al. (2008) in the SDSS mass function 
estimate plotted in Fig. 7. Finally, quite small volumes are 
surveyed when compiling luminosity functions for dwarfs, 
even with the SDSS, so cosmic variance effects may be sig- 
nificant. The discrepancy could then in part reflect differ- 
ences in large-scale structure between the small low-redshift 
volume surveyed by Blanton et al. (2005) and DR7. The 
model also noticeably overpredicts the abundance of very 
luminous galaxies in the (/-band. These are massive galaxies 
undergoing merger-induced starbursts, and it is likely that 
our simple dust model is failing to predict enough extinction 
for these systems. 

4.2 The stellar mass — halo mass relation 

Simplified models for populating dark matter only simula- 
tions with galaxies often assume a simple relation between 
the stellar mass of a galaxy and the mass of the halo sur- 
rounding it - more massive halos should contain more mas- 
sive galaxies at their centres. For such a model to repre- 
sent galaxy clustering even approximately, it must also place 
galaxies at the centres of satellite subhalos, and the resolu- 
tion of the simulation must therefore be good enough that 
a subhalo corresponding to every galaxy can be identified. 
Since tidal stripping often substantially reduces the masses 
of satellite subhalos, but plausibly has little effect on the 
galaxies at their centres, the stellar masses of such galax- 
ies should be much more closely related to the maximum 
masses ever attained by their halos than to their current 
masses. This argument has led many authors to consider 
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Figure 8. Galaxy luminosity functions in the SDSS g, r, i and z photometric bands. The smooth green curves are predictions from our 
preferred model taken from the MS+MS-II at high luminosities and from the MS-II alone at absolute magnitudes fainter than about 
-20. The symbols are observational data for a low-redshift SDSS sample taken from Blanton et al. (2005). 



models which populate simulations with galaxies assuming 
a simple monotonic relation between the stellar mass of a 
galaxy and this maximum past halo mass (Vale & Ostriker 
2004; Kravtsov et al. 2004; Conroy et al. 2006; Wetzel et al. 
2009; Moster et al. 2010; Guo et al. 2010). For cosmological 
simulations of high resolution, matching the (sub)halo abun- 
dance as a function of maximum past mass to the observed 
galaxy abundance as a function of stellar mass allows one 
to derive an (assumed) monotonic relation beween the two 
masses. By using this relation to populate the simulation, 
one can then predict the spatial distribution of galaxies for 
detailed comparison with observation. 

The MS-II simulation provides an unparalleled oppor- 
tunity to carry through this programme because, in combi- 
nation with the MS, it gives a much more precise estimate of 
the abundance of (sub)halos as a function of maximum past 
mass than has previously been available. Guo et al. (2010) 
matched an estimate based on both the MS and the MS-II 
to the SDSS stellar mass function of Li & White (2009), pro- 
ducing the relation between stellar mass and maximum past 
halo mass which we show as a blue curve in Fig. 9. For com- 
parison, green and red symbols show the median value and 
the ±la scatter of stellar mass predicted by our preferred 



model at given past maximum halo mass for z — central 
and satellite galaxies, respectively. Variations in assembly 
history and environmental influence ensure that there is sig- 
nificant scatter in the relation for our model; the rms scatter 
in log M, , is 0.17, 0.20, 0.24 and 0.31 for logM ha i = 14, 13, 
12 and 11 respectively. At low mass, there is a noticeable 
offset between the predictions for satellite and central galax- 
ies, with satellites having systematically larger stellar masses 
for given maximum past halo mass. This behaviour was also 
present in the DLB07 model (see, Wang et al. 2006) and can 
be traced to the fact that low-mass satellite galaxies typi- 
cally achieved their maximum halo mass at z ~ 1, whereas 
for the corresponding centrals this is typically around z ~ 0. 
Since halos are 8 times denser at z = 1 than at z — 0, their 
escape velocities at given mass are roughly 40% larger at the 
higher redshift and this reduces the efficiency with which SN 
feedback can expel gas, increasing the retention of baryons 
for star formation. 

The median stellar mass predicted by our model at each 
maximum past halo mass is very close to the Guo et al. 
(2010) relation at halo masses above about lO 11 A/0, but 
lies noticeably above it at lower masses. This is because Guo 
et al. (2010) extrapolated the Li & White (2009) stellar mass 
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Figure 9. Galaxy stellar mass as a function of maximum past 
halo mass. The latter is the largest mass ever attained by the 
dark matter subhalo centred on the galaxy over its full history. 
This is almost always the mass of the subhalo at the last time its 
central galaxy was type 0, i.e. the present subhalo mass for current 
type galaxies and the subhalo mass just before infall for current 
type 1 and 2 galaxies. Symbols with error bars show predictions 
from our preferred model applied to the MS-II for logAi, < 10. 
and applied to the MS at higher masses. Green symbols are for 
central galaxies (type 0) while red symbols are for satellites (types 
1 and 2). The blue curve is the relation derived directly from the 
SDSS stellar mass function and from subhalo abundances in the 
MS and MS-II under the assumption that the two quantities are 
monotonically related without scatter (Guo et al. 2010). 



function to masses below 1O 83 M using their quoted slope 
of —1.15, which predicts significantly fewer low-mass dwarfs 
than the Baldry et al. (2008) function which we plot in Fig. 7 
and use to set the parameters of our preferred model. In their 
own comparison of a similar relation to observational data on 
satellite galaxy dynamics, More et al. (2009) estimated the 
scatter in log L for relatively massive halos (log M v i r ~ 13) 
to be 0.16 ±0.04. This is in good agreement with the scatter 
actually produced by our galaxy formation model, but is 
considerably smaller than that predicted, for example, by 
the models of Bower et al. (2006) or Font et al. (2008). 



4.3 Gas-phase metal abundances 

In Fig. 10 we show the metallicity of the cold ISM gas as a 
function of stellar mass for star-forming galaxies in our pre- 
ferred model. Here, we define as star-forming those galaxies 
with a specific star formation rate M*/M» > 10~ 11 yr _1 . 
Observational data from Tremonti et al. (2004) and Lee 
et al. (2006) are represented by the solid green curve and 
the red diamonds, respectively. When estimating the oxy- 
gen abundance of model galaxies for comparison with these 
observations, we, for consistency, use the same nucleosyn- 
thetic yields and solar abundances as DLB07. Recent work 
has suggested that these may need to be revised (Asplund 
et al. 2006; Delahaye & Pinsonneault 2006) but the strong- 
line metallicity measurements underlying the observational 
results in Fig. 10 have substantial and controversial uncer- 
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Figure 10. Cold gas metallicity as a function of stellar mass. 
The top panel shows results for star-forming galaxies when our 
preferred model is applied to the MS-II. The bottom panel shows 
similar results but based instead on the DLB07 model. In both 
panels, the solid curves represent observational results for the 
SDSS from Tremonti et al. (2004), while red diamonds are taken 
from Lee et al. (2006) 



tainties, so we prefer to keep our previous assumptions so 
that the models can be easily compared. Results for star- 
forming galaxies in the MS-II are shown as small black dots 
in the upper panel of Fig. 10, which shows that our model 
appears to reproduce the tight observed relation between gas 
metallicity and stellar mass quite well. This is mainly due to 
our introducing a velocity dependence in our SN feedback 
prescription, which leads to less star formation and to more 
effective ejection of metals from low-mass galaxies, thus to 
lower metallicities. 

For comparison, in the bottom panel of Fig. 10 we show 
the predictions obtained when the DLB07 model is applied 
to the MS-II. In this model the SN feedback efficiency is 
assumed to be independent of circular velocity (/3i = /3% = 0) 
leading to a weaker dependence of metallicity on stellar mass 
at low masses than our current preferred model, as well as to 
an overabundance of dwarf galaxies (see Fig. 1). The better 
apparent agreement of the DLB07 model with dwarf galaxy 
properties found in earlier papers turns out to have been 
due largely to the limited resolution of the MS. 
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4.4 Galaxy colours 

The colours of galaxies are influenced strongly by dust, by 
their star-formation histories, particularly by current and 
recent star formation, and by the metallicities of their stars. 
This makes colours especially difficult to predict with models 
of the kind we are discussing, because they are sensitive not 
only to the details of stellar population synthesis, but also 
to assumptions about the production and quantity of dust, 
and about its distribution relative to the different stellar 
populations. While population synthesis models have solid 
theoretical foundations, are well developed and tested, and 
are probably reliable in most situations, the opposite is true 
for dust modelling. For this reason, rather than predicting 
the luminosities and colours of galaxies directly, it is often 
safer to make model predictions for physical properties like 
stellar mass and star formation rate, and to compare these 
with distributions inferred from observation using methods 
designed to be as insensitive as possible to dust. 

In Fig. 11, we show a scatter plot of SDSS u — i colour 
against stellar mass for model galaxies at z = 0. The up- 
per panel includes dust extinction effects, while the middle 
one does not. Blue dots represent galaxies with dominant 
disks (Mbuig e < Mdisk), and red dots galaxies with domi- 
nant bulges (Mbuige > Mdisk). A clear split of the popula- 
tion into a red sequence and a blue cloud is visible in both 
plots. When dust effects are included, our model predicts 
the reddest galaxies to be passive disk systems scattered up 
from the red sequence. It is notable that we predict substan- 
tial numbers of disk galaxies on the red sequence, particu- 
larly at intermediate stellar masses. This appears consistent 
with the fact that SO galaxies substantially outnumber el- 
lipticals in this stellar mass range in the local universe (e.g. 
Dressier 1980), although real SO's rarely have as much dust 
and gas as our model is assigning them. As in DLB07, the 
most massive galaxies are bulge-dominated and lie on the 
red sequence. There are also a few massive bulge-dominated 
galaxies with bluer colours, corresponding to elllipticals that 
have undergone a recent star-formation event, the equivalent 
of the E+A galaxies seen locally (e.g. Zabludoff et al. 1996). 
Finally at low stellar masses we predict both sequences to be 
well populated. As we will see, the fraction of passive dwarf 
galaxies in our model appears larger than observed. To com- 
pare with observation, we show results from SDSS/DR4 in 
the bottom panel. These have been down-sampled to corre- 
spond to a volume-limited subset with stellar masses above 
as in Weinmann et al. (2009). The numbers are 
quite small because the reddest galaxies at this lower mass 
limit fall within the spectroscopic sample only at very low 
redshift. Clear differences with the models appear in the 
slope of the red sequence and in the number of fainter red- 
sequence galaxies. 

To explore this discrepancy with observation in more 
detail, Fig. 12 shows the distributions of u — i (including 
dust extinction) for galaxies in 8 stellar mass ranges span- 
ning four orders of magnitude in stellar mass. The solid 
histograms are constructed from our preferred model ap- 
plied to the MS (for log M t > 10.0) and to the MS-II (at 
lower masses) while the dashed histograms are compiled 
from SDSS/DR7 including 1/V ma x corrections so that they 
correspond to volume-limited statistics. All histograms are 
normalised to have unit integral. For galaxies in the stel- 
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Figure 11. u — i colour as a function of stellar mass for galaxies in 
our preferred model applied to the MS-II. The upper and central 
panels are for model colours including and excluding dust extinc- 
tion effects, respectively. In each panel, red and blue dots refer to 
bulge-dominated and disk-dominated galaxies, respectively, with 
the split set at equal stellar masses for the two components. The 
bottom panel is for a volume-limited subset of SDSS/DR4 with 
no distinction by morphology. 



lar mass range 9.5 < logAf, < 11.0 which contains the 
bulk of all stars, our predictions for the u — i distribution 
are in reasonable agreement with observation, despite our 
over-simplified dust model. At lower masses, the fraction of 
red galaxies is clearly larger in our model than observed. A 
substantial fraction of dwarfs (roughly half) are predicted 
to finish their star formation early and to become passive. 
The observed fraction of such passive dwarfs is substantially 
smaller. At the highest masses, the SDSS galaxies are red- 
der than our model predicts. In the model most of these 
galaxies have mean stellar ages greater than 10 Gyr and stel- 
lar metallicities of order 0.5 Zq. The real galaxies are more 
metal-rich, but for the population synthesis model we are us- 
ing, a 12 Gyr old population with twice solar metallicity has 
u — i = 3.07, thus metallicity and age effects are insufficient 
to explain the discrepancy and no significant dust effects are 
expected. Photometric or K-correction problems may be af- 
fecting these galaxies which are typically at z ~ 0.2. Note 
that at lower mass, the red tails of the distributions corre- 
spond to the (unrealistically) reddened passive disk galaxies 
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Figure 12. u — i colour distributions as a function of stellar 
mass. Solid black curves show the distributions predicted by our 
preferred model (including extinction effects) applied to the MS 
(above logM* = 10.0) and the MS-II (at lower masses), while 
dashed red curves are distributions compiled from SDSS/DR7. 
The range in log M*/Mq corresponding to each panel is indicated 
at top right. 



seen in the upper panel of Fig. 11. This tail is absent at the 
highest masses where the galaxies no longer have gas disks. 



4.5 Tully-Fisher Relation 

There has been a long-standing debate about the ability 
of galaxy formation models in a CDM context to repro- 
duce simultaneously the observed abundance and Tully- 
Fisher (TF) relation of disk galaxies (Kauffmann et al. 1993; 
Cole et al. 1994; Navarro & Steinmetz 2000; Cole et al. 
2000; Blanton et al. 2008). We have shown above that our 
preferred model reproduces the observed galaxy luminosity 
functions in four SDSS bands at z = 0. In this section, we 
study whether it simultaneously produces a relation between 
r-band luminosity and maximum circular velocity which is 
consistent with that observed for isolated spiral and irreg- 
ular systems. Early semi-analytic work on the TF relation 
(e.g. Kauffmann et al. 1993; Cole et al. 1994; Somerville & 
Primack 1999) took the disk rotation velocity to be V v ir, the 
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Figure 13. r-band Tully-Fisher relation. Blue symbols with er- 
ror bars are observational results for isolated disk galaxies taken 
from Blanton et al. (2008) and from Springob et al. (2007). The 
vertical bar on each symbol shows the bin in absolute magnitude 
considered, while the horizontal bar is centred on the median and 
shows the rms scatter of log V ma x for the galaxies within that bin. 
Green dots are results for central (type 0) late-type galaxies from 
our preferred model applied to the MS (brighter than -21) and to 
the MS-II (for fainter galaxies). For the model galaxies Vmax is 
the maximum circular velocity of the hosting dark halo. 



circular velocity of its halo at the virial radius. In the DLB07 
model a reasonable match to the observed relation was in- 
stead found by identifying the disk rotation velocity with the 
maximum circular velocity of its halo (see Croton et al. 2006, 
and its erratum). On the other hand, Cole et al. (2000) found 
that if baryon condensation is assumed to cause halo con- 
traction according to the standard simple formula (Barnes 
1984; Blumenthal et al. 1986) the models are no longer able 
to reproduce the Tully-Fisher relation and the luminosity 
function simultaneously. As discussed in Sec. 3.3, the sim- 
plified model of adiabatic contraction adopted by Cole et al. 
(2000) appears to overestimate the effect of baryons, and 
at least some recent simulations suggest that the maximum 
halo circular velocity found in an equivalent dark matter 
only simulation may be a good approximation to the disk 
rotation velocity (e.g. Tissera et al. 2010). Here we use this 
maximum circular velocity as our disk rotation velocity sur- 
rogate for the TF relation. 

We concentrate on central galaxies in the model and 
compare to observations of isolated systems, because, as 
noted by Blanton et al. (2008) and others Einasto et al. 
(1974), dwarf satellite galaxies appear systematically gas- 
poor and to have systematically lower rotation velocities 
relative to isolated dwarfs of similar stellar mass. This is 
presumably related to the various stripping mechanisms dis- 
cussed above. In order to keep the test simple, it seems wise 
to concentrate on galaxies where such effects are absent. 

We select central galaxies in our model for which the 
r-band absolute magnitude of the bulge is at least 1.5 mag- 
nitudes fainter than that of the galaxy as a whole, and, as be- 
fore, we assume the rotation velocity of the disk to be Knax, 
the maximum circular velocity of its host halo. In massive 
spirals, where baryons dominate in the visible regions, this 
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may underestimate the rotation velocity because we do not 
take the mass of the baryons into account. In dwarf galax- 
ies it may, in contrast, overestimate the rotation velocity 
because baryonic effects are weaker and the observable HI 
may not extend out to the maximum of the halo circular 
velocity curve. For simplicity, we neglect such effects here. 
The Tully-Fisher relation predicted in our preferred model 
by these assumptions is shown using green dots in Fig. 13. 
At absolute magnitudes above -21 the data are taken from 
the MS, while for fainter galaxies they are taken from the 
MS-II. Observational data for relatively bright galaxies from 
Springob et al. (2007) and for isolated dwarfs from Blan- 
ton et al. (2008) are shown by blue symbols. The vertical 
bar on each symbol represents the absolute magnitude bin 
considered and is positioned at the median log Kn ax of the 
observed galaxies in that bin. The horizontal bar shows the 
±1(7 scatter in log V max within the bin. 

It is striking that our model, although clearly not a 
power law, nevertheless agrees reasonably well with the data 
over an absolute magnitude range of about eight magni- 
tudes. There is no evidence for any major problem, even 
for dwarf galaxies with M r ~ —15. This is somewhat un- 
expected, and is due in part to the fact that Blanton et al. 
(2008) excluded dwarf satellite (as opposed to central) galax- 
ies for which the measured rotation velocities are signifi- 
cantly lower at the faintest magnitudes. A more careful com- 
parison does show some discrepancies, however. At high cir- 
cular velocities (V max ~ 250km/s or more) model galaxies 
have a larger scatter in luminosity than the observations. 
The brightest real galaxies have smaller rotation velocities 
than we predict, perhaps because we are stopping star for- 
mation too efficiently in at least some massive systems. 
At the lowest luminosities the simulation predicts slightly 
higher rotation velocities and considerably less scatter than 
is observed. This may reflect the fact that HI data often do 
not reach the peak of the rotation curve in these systems, 
although the current sample of isolated dwarfs is probably 
too sparse to draw reliable conclusions. 

4.6 Profiles and mass functions in rich clusters 

An important aspect of our galaxy formation models is as- 
sociated with the disruption and merging of substructures. 
When tidal effects destroy a dark matter subhalo, we con- 
tinue to follow the properties of its central galaxy, tracking 
its position and velocity using those of the particle which 
was most bound to the subhalo when it was last seen. Such 
"orphan" galaxies may merge with another galaxy (usually 
the central galaxy of the main system) or may themselves 
be tidally destroyed, when specific conditions are satisfied 
(see sections 3.6.2 and 3.7). These procedures account for 
the fact that dark matter subhalos are often prematurely 
disrupted in our simulations both for numerical reasons (res- 
olution may be insufficient to follow tidal stripping down to 
the scale of the central galaxy) and for astrophysical reasons 
(dissipation associated with galaxy formation may make the 
stellar components more resistant to disruption). Thus at 
any given time our galaxy catalogues contain a population 
of orphan (or type 2) galaxies which are concentrated in the 
inner regions of massive halos. 

The large size of our two simulations and the factor 
of 125 difference in their mass resolution makes it possible 
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Figure 14. Projected galaxy number density profiles for samples 
of massive clusters from the MS (red lines) the MS-II (black lines) 
and from the SDSS (blue symbols with error bars). Observational 
and model clusters are selected in the same way and are not scaled 
before stacking (see text for details). Solid lines are for all model 
galaxies with M* > 1.2 X 10 10 Mq, while dashed and dotted lines 
split them into galaxies with surviving dark matter subhalos and 
orphans, respectively. Note the excellent agreement in mean pro- 
file between the MS and MS-II despite the very different number 
of orphans in the two simulations. The SDSS profiles here have 
been corrected for the spectroscopic incompleteness of the sur- 
vey, which varies as a function of projected radius and reaches 
60% near cluster centre. The error bars reflect the uncertainty in 
the mean estimated from the scatter among the 31 SDSS cluster 
profiles. 



to carry out convincing tests of these procedures for the 
first time. Appendix A presents the fraction of galaxies of 
different types in our preferred model. For stellar masses 
in the range 9.5 < log M,/Mq < 11 where both simula- 
tions have good statistics, they show similar fractions of all 
galaxies to be satellites, but the fraction of these satellites 
which are orphans changes from 52% in the MS to 25% in 
the MS-II (at log M»/M = 9.5) or from 27% to 17% (at 
log M*/M@ — 11). Here we test for numerical convergence 
in a considerably more extreme situation by comparing the 
number density profiles predicted for rich clusters in the MS 
and the MS-II. Another sensitive test, based on counts of 
close pairs, is presented below in section 4.9. 

In order to facilitate comparison with real clusters from 
the SDSS, we have implemented a simple "observational" 
cluster finder on our simulations, designed to find objects 
with virial masses in the range 14 < logMvir/M© < 14.5. 
We take all galaxies with stellar masses above 1.2 x 1O 1O M0 
and we view their distribution in "redshift space" where the 
x and y coordinate directions are considered transverse to 
the "line-of-sight" and the z peculiar velocity is added to 
the Hubble constant times the z-coordinate to produce a 
pseudo-recession velocity. We then consider all galaxies as 
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potential cluster centres, and we count neighbours within 
a surrounding cylinder of radius r v = 1.5 Mpc and line- 
of-sight velocity difference ±1200 km/s, weighting by an 
"optimal" filter F(r p ) which we take to be an NFW ap- 
proximation to the projected mass distribution of the target 
clusters. Potential centres are ranked by this weighted neigh- 
bour count and those lying within the cylinder of a higher 
ranked neighbour are eliminated. The MS is then used to 
relate the corresponding unweighted counts N c to halo mass 
in order to identify the count range 45 ^ N c ^ 105 corre- 
sponding to 14 < logM v i r /M0 < 14.5. This algorithm can 
be used almost unmodified on a stellar-mass-limited sample 
of 39 600 SDSS galaxies from DR7 with 0.01 < z < 0.06 and 
M* > 1.2 x 10 10 Af©. The only complication is that the SDSS 
spectroscopy becomes significantly incomplete in the inner 
regions of clusters so that a completeness correction must 
be applied. This can be estimated from the overall spectro- 
scopic completeness as a function of r p within the stacked 
regions. These procedures select 2251, 61 and 31 clusters 
in the MS, the MS-II and the SDSS, respectively 3 . The ef- 
fective SDSS volume searched is 6 x 10 6 Mpc 3 ; given the 
expected cosmic variance expected for the cluster count in 
a volume of this size (~ 25% rms), and the rather large am- 
plitude as = 0.9 adopted in the simulations, the observed 
and simulated cluster abundances appear quite consistent. 

In Fig. 14 we show mean projected number density pro- 
files for stacks of the clusters in these different sets. Solid 
lines show the mean profiles for the two simulations, while 
dashed and dotted profiles split these profiles into galax- 
ies with and without associated dark matter subhalos. Red 
curves refer to the MS and black curves to the MS-II. The 
agreement in the total profiles is remarkable - certainly bet- 
ter than one might have expected since the dashed and dot- 
ted profiles show that orphans make a much larger contribu- 
tion to the MS profiles (where they dominate for r p < 350 
kpc) than to the MS-II profiles (where they dominate only 
for r p < 80 kpc). Within a projected radius of 1.5 Mpc, 
37% of all cluster galaxies more massive than 1O 1O M0 are 
orphans in the MS but only 14% in the MS-II. The fact that 
the total profiles agree so well thus demonstrates that the 
survival times and positions that we assign to our orphans 
are appropriate. 

The SDSS clusters in Fig. 14 are shown by the blue 
symbols with error bars indicating the uncertainty in the 
mean profile due to cluster to cluster variations. The agree- 
ment with the simulations is quite good, although there may 
be an indication that the SDSS clusters are somewhat less 
concentrated than our models. This may be an indication 
that the as value adopted in the simulations is somewhat 
too high (see also Section 4.9 below). 

The biggest halo in the MS-II has a mass of ~ 1O 14 8 M , 
similar to that of the Coma cluster, and contains over 119 
million particles. Its substructures are thus very well re- 
solved. Here we use this biggest halo to investigate whether 
the galaxy stellar mass function inside clusters is expected 
to differ significantly from that of the galaxy population as 
a whole. It is well known that the most massive galaxies 



3 In order to improve the statistics, we include three orthogonal 
projections of the MS-II data, so the mean number of clusters per 
MS-II volume is 20.3. 
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Figure 15. The stellar mass function of galaxies in a rich cluster. 
The solid curve links counts in 0.25dex bins for galaxies within 
_R v ; r = 2Mpc of the the centre of the most massive cluster in the 
MS-II according to our preferred galaxy formation model. Er- 
ror bars indicate Poisson uncertainties in these counts. Red open 
triangles represent the general stellar mass function of galaxies 
constructed from the MS-II as a whole. This has been renormal- 
ized arbitrarily to allow its shape to be compared to that of the 
cluster stellar mass function. 



occur exclusively in rich clusters, and that cluster popula- 
tions have systematically different star formation histories 
and morphologies to field galaxies. Nearby clusters also ap- 
pear to contain a population of small dwarf ellipticals which 
are not found in less dense environments (e.g. Binggeli et al. 
1990). Thus it is interesting to see whether our galaxy for- 
mation model predicts differences which might correspond 
to these observations, and, in particular, to see if the rela- 
tive number of dwarf galaxies in a rich cluster is predicted 
to differ from that in the "field". 

We study this in Fig. 15. The solid curve is the stellar 
mass function for galaxies within R V [ T = 2Mpc of the centre 
of this massive cluster, with error bars indicating the Poisson 
uncertainty in the count in each bin. The slope at the low 
mass end is around -1.4, which is higher than the observed 
r- or _R-band slope for galaxies in the Coma cluster: ~ 1.16 
(Beijersbergen et al. 2002; Mobasher et al. 2003), but per- 
haps consistent with recent observational estimates based on 
the SDSS data for nearby X-ray-selected clusters (Popesso 
et al. 2006). At very faint magnitudes the slope in the Coma 
cluster may be steeper (Adami et al. 2007; Jenkins et al. 
2007; Milne et al. 2007). Given the large dispersion in ob- 
servational results, our model seems quite compatible with 
the data. The triangles show the overall stellar mass function 
of the MS-II, renormalized for ease of comparison with the 
cluster result. The shapes of the two stellar mass functions 
are very similar, both the faint-end slope and the break at 
high mass. This echoes the results found for the infrared lu- 
minosity function of the Coma cluster by Bai et al. (2006). 
This is interesting, since both observations and the simu- 
lations of this paper show substantial differences in colour 
and morphology between clusters and the field. In the simu- 
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Figure 16. The stellar mass fraction in intergalactic stars as a 
function of virial mass for clusters. The solid black line shows 
the fraction of all stars within R v \ r which are assigned to the 
intracluster component when our preferred model is applied to 
the MS. Dashed black lines show the 16 and 84% points of the 
distribution of this fraction. Solid and dashed red lines show the 
same statistics but for the fraction of stars in the main subhalo 
of each cluster which are associated with its diffuse component, 
rather than with its central galaxy. 



lations over 95% of cluster galaxies within R vil are passive. 
This fraction seems overly large in comparison to observa- 
tion (e.g. Hansen et al. 2009), again reflecting the fact that 
the passive galaxy fraction in general is somewhat too high 
in our model. 



4.7 Intracluster Light 

Recent observations of diffuse intracluster light and of intra- 
cluster stars (Zibetti et al. 2005; Gerhard et al. 2005; Mihos 
et al. 2005; Gonzalez et al. 2005; Aguerri et al. 2006; Gonza- 
lez et al. 2007; McGee & Balogh 2010) indicate that a signif- 
icant fraction of all cluster stars lie between the galaxies, but 
they disagree about the exact amount. It seems likely that 
such stars must be the remains of disrupted galaxies, and our 
model now includes a treatment of the tidal disruption pro- 
cess. In Fig. 16, we show the fraction of cluster stars in the 
intergalactic component as a function of cluster virial mass. 
We consider two different fractions here. The black lines re- 
fer to the fraction by mass of all stars within 7? v ir which are 
assigned to the intergalactic component. The solid curve is 
the median value at each M v i r , while the dashed lines indi- 
cate the 16 and 84% points of the distribution. This intra- 
cluster fraction increases with cluster mass and has a large 
scatter in low-mass clusters. In our preferred model (here 
applied to the MS) around 5-10% of all stars in clusters 
with M v i r > 5 x 1O 14 M ) are in the intracluster component 
and the dependence on cluster mass is quite weak. In less 
massive systems this fraction drops very rapidly, reaching 
1% in groups of mass 3xlO 13 M0. Both the trends and the 
value are within the scatter of the observational results cited 
above. 

Fig. 16 also shows another fraction of interest. The red 



curves show the median and the 20 and 80% points of the 
distribution of the fraction of all the stars in the main sub- 
halo which are associated with the diffuse component, rather 
than with the central galaxy. This can be considered as a 
proxy for the fraction of the stellar mass of the cD galaxy 
which is associated with its extended envelope. This frac- 
tion also increases with cluster mass, ranging from ~ 10% 
in clusters with M v i r ~ 1O 14 M0 to 30% in clusters with 
M v i r ~ 1.4 x 10 M . Thus, in the richest clusters, the mass 
in intergalactic stars is comparable to the stellar mass of 
the main body of the central galaxy, or, alternatively, the 
extended envelope of the cD galaxy contains about half of 
its stars. In galaxy groups, this fraction decreases rapidly 
with decreasing virial mass, reaching 1% in groups of mass 
2xlO 13 M . 

4.8 Luminosity function of Milky Way satellites 

The abundance of the very lowest mass galaxies can be mea- 
sured observationally only in the Local Group, in particular, 
in the halo of the Milky Way. The apparent discrepancy be- 
tween the relatively small number of observed satellites and 
the many dark matter subhalos predicted in a ACDM cos- 
mogony has been promoted as "the missing satellite prob- 
lem" , a possible flaw in the concordance structure formation 
model (Moore et al. 1999; Klypin et al. 1999), despite earlier 
suggestions that it might rather reflect the astrophysics of 
galaxy formation in weak potential wells (Kauffmann et al. 
1993). Over the last decade new observational results, pri- 
marily from the SDSS, have increased the directly observed 
number of satellites by almost a factor of two and the esti- 
mated total number of satellites by about a factor of four 
(e.g. Koposov et al. 2008). At the same time, improved sim- 
ulations have increased the predicted number of subhalos by 
a factor of 1000 (e.g. Springel et al. 2008) . Thus the discrep- 
ancy has grown. Our galaxy formation models make it pos- 
sible to address this issue in the context of the more general 
problem of matching the low-mass end of the stellar mass 
function of galaxies. This is because the MS-II contains sev- 
eral thousand isolated galaxies similar in mass to the Milky 
Way, and its resolution turns out to be (just) sufficient to 
get predictions for objects with stellar masses comparable 
to those of the observed Milky Way stellites. 

In the MS-II, there are around 7000 halos with virial 
mass within a factor of three of that estimated for the halo 
of the Milky Way (see Boylan-Kolchin et al. (2010) for an 
analysis of the properties of these halos and their substruc- 
ture). In order to make a detailed comparison, we select all 
disk-dominated (M*,disk > M*, bulge) central galaxies with 
total stellar mass between 4 and 8 x10 10 Mq. (The stellar 
mass of the Milky Way is estimated to be 5 x 1O 1O M (Flynn 
et al. 2006).) This provides us with a sample of 1603 "Milky 
Ways" which have median halo mass M vir = 1.30x 1O 12 M 
with lower and upper quartiles at 0.90 and 2.18 xlO 12 M . 
For the purposes of this section, all galaxies within 280 kpc 
of each "Milky Way" are defined to be its satellites. Fig. 17 
shows the cumulative ^-band luminosity function of these 
satellite systems in our preferred model and in two varia- 
tions with different assumptions about reionization. Specif- 
ically, we plot the median and the 10 and 90% points of 
the distribution of satellite counts as a function of limiting 
absolute magnitude, My. A dashed red curve plotted for 
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Figure 17. Cumulative luminosity functions for the Milky Way 
satellite system, defined to consist of all galaxies within 280 kpc 
of the Galactic Centre. Simulated "Milky Ways" are taken to 
be disk-dominated central galaxies with stellar masses between 
4 and 8x10 10 Mq. Solid curves give the median satellite count 
predicted above each absolute magnitude, while dotted curves 
delineate the 10% and 90% tails of the count distribution. The 
upper panel gives results for our preferred model applied to the 
MS-II. This assumes the effects of reionization to be as advo- 
cated by Okamoto et al. (2008). In the central panel we show 
what happens if we instead use the reionization prescription of 
Gnedin (2000), keeping all other model parameters fixed. Reion- 
ization effects are weaker in the Okamoto et al. (2008) model than 
in that of Gnedin (2000). More detailed discussion of these two 
recipes can be found in Sec. 3.1. For the lower panel, reionization 
is assumed to have no effect on galaxy formation. In each panel 
the cumulative luminosity function for the 11 "classical" satel- 
lites of the Milky Way is shown as a stepped red curve ending at 
My ~ — 8. The abundance of satellites with My < — 5 estimated 
by Koposov et al. (2008) is indicated by a large filled red circle. 
Because of the substantial and uncertain completeness correction 
needed to make this estimate, we have arbitrarily assigned it an 
error bar of a factor of two. 
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Figure 18. The effects of reionization on the low-mass end of the 
stellar mass function of galaxies. The red curve is the ratio of the 
stellar mass function predicted for the MS-II by a model exclud- 
ing the effects of reionization to that predicted by our preferred 
model which is identical except that reionization is included fol- 
lowing the prescription of Okamoto et al. (2008). Reionization 
changes the abundance of galaxies only at stellar masses below 
10 8 Mq. Effects are stronger if the prescriptions of Gnedin (2000) 
are used instead, as in DLB07. This is shown by the green curve 
which gives the ratio of the abundances predicted for this model 
to those predicted by our preferred model. Above 10 8 Mq the 
effects remain below 20%. 



Mv < — 8 represents the cumulative luminosity function of 
the 11 "classical" Milky Way satellites. To this limit, the ob- 
served sample is thought to be (almost) complete. We also 
use a large filled red circle to indicate the estimate of 45 
Milky Way satellites with My < -5 and r < 280 kpc from 
Koposov et al. (2008). This estimate required a large and 
uncertain incompleteness correction, so we have arbitrarily 
assigned it an error bar of a factor of two. 

The top panel of Fig. 17 shows results for our pre- 
ferred model which assumes the Okamoto et al. (2008) pre- 
scriptions when estimating the effects of reionization. The 
predicted satellite abundance is consistent with observation 
all the way from bright LMC/M33-like systems down to 
Mv ~ —11, even though model parameters were set to 
match the general galaxy stellar mass function rather than 
Local Group data. For fainter systems, the observational es- 
timates are close to the lower 10% point of the predicted 
counts, but, as just noted, the Koposov estimate has a sub- 
stantial intrinsic uncertainty. In addition the classical satel- 
lite count may well have missed a one or two systems be- 
hind the Galactic Plane. As the middle panel shows, if we 
substitute the Gnedin (2000) parameters used by DLB07 
for those of Okamoto et al. (2008), the predicted number 
of faint galaxies is reduced, and the match to the observa- 
tional estimates is almost perfect. However, Okamoto et al. 
(2008) and Hoeft et al. (2006) argue that the simulations 
of Gnedin (2000) substantially overestimated the extent to 
which an ionizing background suppresses the accretion of gas 
onto small haloes. If, on the other hand, we neglect the ef- 
fects of reionization altogether, the bottom panel shows the 
disagreement with the observational data to worsen only at 
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the faintest magnitudes. The median count of satellites with 
My < —5 is predicted to be about four times the Koposov 
estimate, but brighter than My ~ —10, the abundances 
are almost unchanged from our preferred model. Thus, if 
Okamoto et al. (2008) are right, reionization has a significant 
effect only on the very faintest galaxies. This is consistent 
with results of previous work (Bullock et al. 2000; Somerville 
2002; Benson et al. 2003; Gnedin & Kravtsov 2006; Okamoto^ 
et al. 2010) iL 
We explore this point further in Fig. 18, which shows £ 
how reionization modelling affects the low-mass end of the * 
overall stellar mass function. We plot the factor by which the 
galaxy abundance in the MS-II is changed as a function of 
stellar mass if our preferred model, which uses the Okamoto 
et al. (2008) reionization parameters, is altered to use those 
of Gnedin (2000), as in DLB07 (green line), or to neglect 
the effects of reionization altogether (red line). In our pre- 
ferred model, reionization affects the abundance of galaxies 
noticeably only below about 10 7 A/q. The stronger effects 
implied by the Gnedin (2000) recipe, reduce the abundance 
by about 20% already at 1O 8 M0, but remain small for more 
massive systems. Thus we conclude that reionization has very 
little effect on galaxies similar to the brighter Local Group 
dwarfs, but may significantly affect the abundance of the 
fainter dwarf spheroidals. 

4.9 Correlation Functions 

The SDSS has revolutionised our knowledge of the nearby 
galaxy population not only by providing quantitatively reli- 
able data for galaxy abundances as a function of luminosity, 
stellar mass and colour over the full range from dwarfs to 
cD galaxies, but also by providing precise measurements of 
the spatial clustering of galaxies as a function of their lu- 
minosity and colour on scales from 20 kpc to 30 Mpc and 
beyond. With simulations the size of the MS and the MS-II, 
our galaxy formation models make equally precise predic- 
tions for the clustering of simulated galaxies as a function of 
their physical properties. Comparing observation and simu- 
lation provides powerful constraints on the galaxy formation 
modelling. No modern semi-analytic or hydrodynamic sim- 
ulation of the formation of the galaxy population should be 
considered viable unless it demonstrates at least adequate 
agreement, not only with stellar mass, luminosity and color 
distributions, but also with clustering as a function of galaxy 
properties. 

In Fig. 19 we compare the projected autocorrelation 
of stellar mass in the final release of the SDSS to the re- 
sults we obtain for our preferred galaxy formation model. 
In the upper panel, red and blue symbols are results from 
the MS and MS-II, respectively; the black solid line shows 
the SDSS/DR7 measurement from Li & White (2009). The 
error bars on the latter include the effects of counting noise 
and cosmic variance and are impressively small. This is re- 
markable because small-scale correlations are dominated by 
the distribution of satellite galaxies near halo centre, where 
one might have expected resolution effects to cause substan- 
tial differences. For example, the number of type 2 (orphan) 
galaxies differs substantially between the two simulations 
(see Appendix A). In part, the agreement reflects the fact 
that, as Li & White (2009) show, the main contribution to 
the autocorrelation of stellar mass comes from galaxies with 



1000 



1 1 1 1,1 

r- » . 




: t • • 




i i ... 


... i . . . L . . 


^ — . — 




-= 



0.01 0.10 1.00 10.00 100.00 

r p [Mpc] 



Figure 19. The projected autocorrelation function of stellar mass 
(upper panel). Blue and red circles show results from our preferred 
model applied to the MS-II and to the MS respectively. Numerical 
convergence is excellent, even on scales below 100 kpc. An esti- 
mate from the final release of the SDSS is shown by a black solid 
line joining points with error bars which include both counting 
noise and cosmic variance (Li & White 2009). On large scales our 
model overstimates the observed amplitude of clustering by 10 to 
20%. On small scales the discrepancy rises to a factor of two. In 
the lower panel we show the ratio of the two model autocorrela- 
tion functions to the SDSS estimate. 



individual stellar masses similar to the Milky Way, and thus 
well above the resolution limit of the MS (see the stellar 
mass functions in Sec. 4.1 and the mass-dependent corre- 
lation functions presented below). For r p > 2 Mpc, where 
the correlations are produced by galaxies inhabiting differ- 
ent halos (thus typically both type galaxies), the model 
autocorrelation function is 10 to 20% higher than that ob- 
served. The small difference between the MS and the MS-II 
on these scales may reflect the effect of cosmic variance due 
to the relative small box size of the MS-II. On smaller scales 
where the correlations are dominated by galaxy pairs inhab- 
iting the same halo (thus typically type - type 1, or type 
- type 2 pairs) the discrepancy grows, reaching a factor of 2 
at r v < 100 kpc. This suggests an overdominance of 1-halo 
relative to 2-halo pairs in comparison to the observations, 
arguing, perhaps, for a lower value of as than used in the 
MS cosmology (see Li & White 2009). 

We investigate the source of this discrepancy further in 
Fig. 20, which shows projected autocorrelation functions for 
galaxies in a set of disjoint stellar mass ranges, as indicated 
by the labels in each panel. Black solid and blue dashed 
curves give the predictions obtained by applying our pre- 
ferred galaxy formation model to the MS and to the MS-II, 
respectively. Corresponding 1/V max - weighted estimates from 
the full SDSS /DR7, obtained using the techniques of Li et al. 
(2006), are shown by symbols with error bars. Here the er- 
rors are estimated from a set of 80 mock SDSS surveys and 
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Figure 20. Projected autocorrelation functions for galaxies in different stellar mass ranges. Black solid and blue dashed curves give 
results for our preferred model applied to the MS and the MS-II, respectively. Symbols with error bars are results for SDSS/DR7 
calculated using the same techniques as in Li et al. (2006). The two simulations give convergent results for Af* > 6 X 10 9 Mq. At lower 
mass the MS underestimates the correlations on small scales but still matches the MS-II for r p > 1 Mpc. The model agrees quite well 
with the SDSS at all separations for M* > 6 X 10 10 Mq, overestimating the correlations slightly on small scales, but at smaller masses 
the correlations are overestimated substantially, particularly at small separations. 



so should, in principle, include cosmic variance effects. This 
becomes a significant issue at the smallest masses. No re- 
sult is shown for the MS-II in the most massive bin, because 
it contains too few galaxies to give a meaningful estimate. 
Results from the two simulations converge for galaxies more 
massive than 6 x 1O 9 M0. For smaller masses the MS un- 
derpredicts the correlations on small scales but still agrees 
with the MS-II for r p > 1 Mpc. This indicates that resolu- 
tion limitations begin to affect satellite galaxies in the MS 
at higher stellar mass than central galaxies. 

For M* ^ 6x 1O 1O M0 the model autocorrelations agree 
with the SDSS at all separations to better than about 20%. 
For M* > 6 x 1O 9 M0, simulation and observation con- 
tinue to agree at about the 20% level for r p > 2 Mpc. 
This shows that the relation between halo mass and cen- 
tral galaxy mass shown in Fig. 9 leads to autocorrelations 
for central galaxies as a function of their stellar mass which 
are in good agreement with observation. The small remain- 
ing off-set may indicate a fluctuation amplitude somewhat 
smaller than the as = 0.9 adopted in the simulations. At 
yet smaller masses the large-scale correlation amplitude es- 
timated from the SDSS disagrees with the model. Plots of 



the distribution of these galaxies on the sky show that their 
correlations are dominated by a very small number of struc- 
tures (just the Coma and Virgo clusters in the lowest mass 
bin) which are particularly pronounced in the minority red 
population. In these very shallow samples, correlation es- 
timates are also significantly distorted by peculiar velocity 
effects (e.g. the finger-of-god of the Coma cluster and Virgo- 
centric infall). Proper accounting for these effects is beyond 
the scope of this paper. 

At smaller separations (r p ^ 1 Mpc) Fig. 20 shows 
substantial discrepancies between model and observation for 
stellar masses below 6 x 10 10 Mq, indicating that there are 
more satellite-central pairs in the model than in the real 
data. Since the overall abundance of galaxies as a function 
of stellar mass matches observation very well (see Fig. 7), 
this discrepancy indicates that too large a fraction of the 
model galaxies are satellites. Again this is a clear indication 
favoring a lower value of as which would result in a lower 
abundance of the high-mass halos which host two or more 
galaxies in these stellar mass ranges (cf van den Bosch et al. 
2007). 

Additional insight into possible errors in our treatment 
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Figure 21. Projected autocorrelation functions for galaxies as a function of colour and stellar mass. As in Fig. 20, solid and dashed curves 
are for our preferred model applied to the MS and to the MS-II, respectively. Symbols with error bars are again derived from SDSS/DR7 
using the techniques of Li et al. (2006). In each mass range, the galaxies are split into passive (red) and active (blue) subsamples according 
to their g — r colour. The colours of the symbols and curves correspond to those of the populations. Qualitatively, the agreement between 
models and observations is good, with quantitative agreement at both high (M* > 6 X 1O 1O M0) and low (M* < 6 X 10 9 Mq) stellar mass 
and a somewhat stronger dependence of clustering on colour than is observed at intermediate stellar masses. 



of the astrophysics of galaxy evolution can be obtained by 
studying clustering as a function of star formation activity. 
To this end, Fig. 21 repeats Fig. 20 but with the galaxies 
in each mass range divided into "passive" (red) and "ac- 
tively star-forming" (blue) subsamples according to their 
g — r colours, as in Li et al. (2006). 4 Lines and symbols 
are as in Fig. 20, except that they are coloured according 
to the colour of the corresponding galaxy population. As 
expected, red galaxies are more clustered than blue galax- 
ies on all scales and at all stellar masses. It is encouraging 
that the effects are qualitatively similar in the models and 
in the SDSS data. Indeed, at large separation (r p > 2 Mpc) 
there is reasonable quantitative agreement for both popu- 
lations at all but the smallest stellar masses, while at large 
stellar mass (M« ^ 6 x 1O 1O M ) there is good agreement 
at all separations. For active galaxies, this simply indicates 
once more that our halo mass - central galaxy mass rela- 
tion leads to the right large-scale correlations as a function 



4 For the simulations we take the division at the minimum of the 
"green valley" in a plot similar to Fig. 11. 



of M* for type galaxies. For passive galaxies the situa- 
tion is more complex, since most of the lower mass objects 
are satellites rather than centrals. Apparently, at given stel- 
lar mass, their distribution across halos of different mass is 
similar in the simulation and in the real world. For the two 
lowest mass bins the large-scale correlations are again dis- 
torted by the small volume and peculiar velocity distortion 
effects discussed above. 

At small separations the simulations overpredict the au- 
tocorrelations of passive galaxies for stellar masses in the 
range 6 x 1O 9 M < M, < 6 x 10 10 A/ Q , but, curiously, the 
MS-II again matches the real data at lower mass. Small- 
scale correlations of active galaxies are underpredicted in 
all our lower stellar mass bins, showing that our model still 
has somewhat too few blue satellite galaxies. An interest- 
ing example is provided by our lowest stellar mass bin. For 
2 Mpc > r p > 200 kpc the MS-II model fits the SDSS data 
quite well in Fig. 21, yet lies substantially above them in 
Fig. 20. This is because our model overpredicts the fraction 
of passive galaxies in this mass range (see Fig. 11). Similar 
apparent discrepancies between the model/observation com- 
parison in Fig. 21 and that in Fig. 20 are visble at a lower 
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level in other mass ranges, and again reflect the slightly dif- 
ferent weightings applied in the two cases when going from 
colour-differentiated to "total" results. 

4.10 Some properties at higher redshift 

So far we have only discussed properties of our models at 
z ~ 0. This is because our observational knowledge of the 
galaxy population is still much more complete, more precise 
and less subject to systematic error in the nearby universe 
than at high redshift, despite the enormous recent progress 
in amassing data for relatively large, objectively selected 
samples of distant galaxies. Nevertheless, a viable galaxy for- 
mation model must be consistent not only with the present- 
day galaxy population, but also with that at all earlier times, 
so a comparison of our models with high-redshift popula- 
tions is a critical part of assessing how realistically they treat 
the astrophysics of galaxy formation. Such work is compli- 
cated by the strong selection effects and the substantial ob- 
servational uncertainties which affect the measurement of 
physical properties for faint and distant galaxies. As a re- 
sult, detailed comparison is beyond the scope of the present 
paper. Earlier versions of our models have been compared 
to the evolution of the cosmic star-formation rate density 
by Croton et al. (2006), to the evolution of brightest cluster 
galaxies out to z = 1 by De Lucia & Blaizot (2007), to the 
galaxy counts, luminosity functions and redshift distribu- 
tions inferred from deep magnitude-limited redshift surveys 
by Kitzbichler & White (2007) and to the abundances, red- 
shift distributions, stellar mass distributions and clustering 
of colour-selected samples of z ~ 2 and z ~ 3 galaxies by 
Guo & White (2009). The current models can be expected to 
give similar results to this previous work and to be sensitive 
to many of the same uncertainties, notably to the treatment 
of dust obscuration. In this section we will limit ourselves to 
presenting two of the least uncertain model predictions at 
high redshift. 

In Fig. 22 we compare the evolution of the cosmic star 
formation rate density predicted by our preferred model to 
a compilation of observational estimates taken from Hop- 
kins (2007). The most obvious feature of this plot is a clear 
off-set between the model and the observations. At all red- 
shifts the model lies a factor of two or more below the centre 
of the cloud of observational points. This is a reflection of 
the well known fact that if one integrates observational esti- 
mates of the star formation rate density with respect to time, 
one substantially overpredicts the observed stellar mass den- 
sity, not only at z — but also at all higher redshifts (e.g. 
Wilkins et al. 2008). We have chosen to adjust our model to 
fit the SDSS stellar mass function, so we necessarily fail to 
fit observational estimates of the evolution of the star for- 
mation rate density. In our model the rate of star formation 
peaks at z ~ 3 and has already declined again by a fac- 
tor of 3 at z ~ 1, whereas the observations suggest a more 
constant star formation rate density over this time interval. 
Given the large scatter in the observational estimates and 
the discrepancy just discussed, it is difficult to know how 
seriously to take this difference. As we shall see in the next 
paragraph, however, there are other indications that galaxy 
formation occurs too early in our model, particularly for 
low-mass galaxies. 

Stellar masses for high-redshift galaxies are notoriously 



difficult to estimate because of the faintness of the images, 
the strong effects of dust, and the fact that the observed op- 
tical and near-IR bands correspond to the rest-frame blue 
and ultraviolet. The situation has improved considerably 
with the availability of deep data at 3.6 to 8/u from Spitzer, 
and according to the careful error analysis of Marchesini 
et al. (2009), masses with realistic error bars can now be 
estimated out to at least z ~ 4. In Fig. 23 we compare 
the stellar mass functions predicted by our preferred model 
to recent observational estimates based on combined very 
deep optical, near-IR and Spitzer photometry from Perez- 
Gonzalez et al. (2008) and Marchesini et al. (2009). We have 
shifted all these observational estimates so that they corre- 
spond to the same Chabrier Initial Mass Function used in 
our models. As Marchesini et al. (2009) describe, even with 
this excellent data coverage substantial random errors re- 
main in the stellar masses estimated for individual galaxies 
(see also Fontanot et al. 2009). To account roughly for this, 
we convolve the stellar mass functions predicted by our pre- 
ferred model with a gaussian of dispersion 0.25 dex in log M* 
before comparing them with the observations. 

Our model parameters are adjusted so that they fit the 
observed stellar mass function at z ~ 0. This good agree- 
ment is maintained out to redshifts somewhat less than 
unity. At higher redshift, the massive end of our predicted 
mass functions remains consistent with observation, once 
it has been convolved with the observational mass estima- 
tion uncertainties, but the abundance of lower mass galax- 
ies ( M* ~ 10 10 Mq) is substantially overpredicted. 5 At face 
value, the discepancy suggests that low-mass galaxies form 
considerably earlier in our model than in the real universe. 
This is consistent both with the overly high redshift of the 
peak of the star formation rate density (see Fig. 22) and the 
overly large fraction of passive galaxies in the z ~ low- 
mass population (see Fig. 12). The problem is not specific to 
the details of our model. It has been seen in the comparison 
of earlier models (both our own and those of others) to this 
and other similar datasets (e.g. Fontana et al. 2006; March- 
esini et al. 2009; Lo Faro et al. 2009; Fontanot et al. 2009; 
Cirasuolo et al. 2010). As several of these authors empha- 
sise in their own discussion, the problem seems most likely 
to lie in the way star formation is treated in the models, 
particularly at high redshift. 



5 DISCUSSION 

New observational data at low redshift give precise measures 
of the abundance and clustering of galaxies as a function of 
their physical properties (stellar mass, luminosity, size, star 
formation rate, nuclear activity...) over a range of almost five 
orders of magnitude in stellar mass (7 < log M*/M Q < 12). 
Abundances of even lower mass galaxies are measured rea- 
sonably reliably in the Local Group. In addition, the explo- 
sion of data from ultra-deep surveys is beginning to provide 
convincing results for the general galaxy population at much 

If the systematic error ranges discussed by Marchesini ct al. 
(2009) are considered appropriate, this overprediction appears 
only marginally significant, but a large part of these systematic 
errors are due to possible IMF variations which we exclude for 
the present discussion. 
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Figure 23. Stellar mass functions for a series of redshift intervals indicated by the labels in each panel. Observational data are taken 
from Perez-Gonzalez et al. (2008) and from Marchesini et al. (2009). Marchesini et al. (2009) compiled their mass functions for wider 
bins than Perez-Gonzalez et al. (2008) so in each panel we plot the Marchesini et al. (2009) results for the wider bin that includes 
the indicated redshift range. For the triangles representing the Perez-Gonzalez et al. (2008) data we use the error bars quoted in their 
paper. For the filled circles representing the Marchesini et al. (2009) results we use the error estimates which include counting statistics, 
cosmic variance, photometric redshift uncertainties and photometric errors, but exclude systematic uncertainties due to the IMF and 
other stellar population modelling issues. The mass scales of these observational results have been shifted to correct approximately to the 
Chabrier IMF assumed in our modelling. Black curves are the functions measured directly from the MS and the MS-II for our preferred 
galaxy formation model, while green curves show the result of convolving with a gaussian of dispersion 0.25 dex in logM, in order to 
represent uncertainties in the individual observational stellar mass determinations. 



earlier cosmic epochs. Matching such a wealth of data over 
such a large dynamic range is an extraordinary challenge for 
any a priori galaxy formation model. By combining results 
from the MS and the MS-II, and by updating and readjust- 
ing our treatments of the many relevant astrophysical pro- 
cesses, we have made a model which has the necessary dy- 
namic range and statistical power to confront the full range 
of abundance and clustering data available at low redshift. 
The MS-I gives good statistics for rare, high-mass galaxies, 
while the MS-II provides well-resolved assembly histories for 
low-mass systems. 

In this paper we have extended and modified our ear- 
lier treatments of the transition between the rapid infall 
and cooling flow regimes of gas accretion, of the sizes of 
bulges and of gaseous and stellar disks, of supernova feed- 



back in low-mass galaxies, of the transition between central 
and satellite status as galaxies fall into larger systems, and 
of the stripping of gas and stars once they have become 
satellites. For physically plausible values of its parameters, 
the new model fits both the abundance and the large-scale 
clustering of low-z galaxies as a function of stellar mass, lu- 
minosity and (to a lesser extent) colour. At high mass the 
efficiency of star formation is limited by AGN feedback, as 
proposed by Croton et al. (2006). At low mass, consistency 
with the observed SDSS luminosity and stellar mass func- 
tions requires supernova feedback to be significantly more 
efficient and the reincorporation of ejected gas to be con- 
siderably less efficient than in DLB07. This enhanced SN 
feedback also leads to reasonable agreement with the abun- 
dance of faint satellites around the Milky Way, suggesting 
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Figure 22. Cosmic star formation rate density as a function of 
rcdshift. The crosses are individual observational estimates com- 
piled by Hopkins (2007) while the solid curve is obtained from 
our preferred model applied to the MS. 



that reionisation influences the formation of, at most, the 
very smallest galaxies (see also Li et al. 2010; Maccio et al. 
2010). 

For galaxies of high stellar mass, our preferred model 
also fits both the colour distribution and the small-scale clus- 
tering of SDSS galaxies (logM*/Mo > 9.5 for the colours 
and > 10.5 for the clustering). At lower stellar mass, the 
model predicts a substantial fraction of red, passive galaxies 
which are not present in the SDSS data, and a clustering 
strength which rises progressively above that observed for 
r p < 1 Mpc. (Note, however, that the difference in cluster- 
ing between active and passive galaxies is still modelled quite 
accurately.) Given that our model matches both the stellar 
mass function and the mass-dependent large-scale clustering 
data from SDSS, this excessive small-scale clustering im- 
plies that too large a fraction of our galaxies are satellites 
at each stellar mass. Since individual groups and clusters 
in our model have galaxy occupation numbers and radial 
distributions in quite good agreement with observation, the 
discrepant small-scale correlations suggest that massive ha- 
los are overabundant in our simulations, i.e. that as = 0.9 is 
too large (c.f. van den Bosch et al. 2007). We intend to test 
this explicitly in future work by using the rescaling tech- 
niques of Angulo & White (2010) on the MS and MS-II so 
that they can be used to construct galaxy formation mod- 
els very similar to those of this paper, but for cosmologies 
other than that originally assumed for the simulations, for 
example, cosmologies with lower values of as, as suggested 
by more recent WMAP results. 

The excessive passive fraction at low stellar mass im- 
plies that our preferred model is quenching star formation 
in small halos in order to limit the total production of stars, 
whereas real objects form stars at a steady but low rate un- 
til the present day. This is also the principal reason why 
the model continues to have too few blue satellites, despite 
our improved treatment of stripping effects - at low stellar 
masses (log M,/Mq < 10) there are too few star- forming 
galaxies everywhere. Low-mass star-forming galaxies in the 



model fit on the observed Tully-Fisher relation for isolated 
galaxies just as well as their giant cousins, and their large- 
scale clustering is also correct. Thus dwarfs appear to be 
forming in the proper dark halos. The overly early trun- 
cation of their star formation is very likely related to the 
fact that while the model correctly fits the observed abun- 
dance of massive galaxies (M* ~ 10 11 Mq) out to z ~ 4, it 
overpredicts the observed abundance of lower mass systems 
(M, ~ lO lo Ai0) by progressively larger amounts beyond 
z ~ 0.6. Lower mass galaxies clearly complete their forma- 
tion too early in the model. 

With the increased resolution provided by the MS-II we 
are able to show that the stellar mass function of galaxies in 
rich clusters is predicted to be very similar in shape to that 
in the general field, even down to M* ~ 1O 7 M . Almost 
all galaxies within the virial radius of a relaxed cluster are 
predicted to be passive, but this may be an overestimate 
for the reasons discussed in the last paragraph. Our new 
treatment of galaxy disruption suggests that 5% to 10% of 
all cluster stars should be be part of the intracluster light, 
and that this fraction should increase with cluster mass and 
show substantial cluster to cluster variation. 

The predictions for the luminosity functions and ra- 
dial number count profiles of clusters are very similar in 
the MS and the MS-II and agree quite well with observa- 
tion provided the substantial population of galaxies with- 
out surviving dark matter subhalos is included. Such orphan 
galaxies account for almost half of all cluster members with 
M, > 10 10 M» in the MS, and for about 13% in the MS- 
II. Without them the abundance of galaxies in the inner 
cluster would be substantially underpredicted. This demon- 
strates that, even at MS-II resolution, schemes that place 
galaxies in subhalos in a high-resolution simulation without 
accounting for subhalos which have been tidally disrupted 
but whose galaxies have survived (e.g. Vale & Ostriker 2004; 
Conroy et al. 2006; Wetzel et al. 2009; Moster et al. 2010; 
Guo et al. 2010) will not correctly reproduce the observed 
structure of galaxy clusters. This argument was already pre- 
sented by Gao et al. (2004). 

The degree to which our physically based model repro- 
duces the observed abundance and clustering properties of 
the z ~ galaxy population is impressive, but there are 
clear and significant discrepancies, and a comparison with 
high-redshift populations, although barely started here, also 
shows substantial discrepancies. Further work is needed to 
understand the source of these problems. Our simple recipes 
for complex astrophysical processes may turn out to be in- 
appropriate when more accurate treatments become feasi- 
ble. In addition, processes other than those we discuss may 
produce similar behaviour, making them operationally in- 
distinguishable at the present level of description. Finally, 
there are undoubtedly degeneracies among the model pa- 
rameters we have adjusted, making our specific model non- 
unique (see, for example, Henriques et al. 2009; Bower et al. 
2010; Neistein & Weinmann 2009). Such degeneracies can 
only be lifted, and the recipes improved, by increasing the 
range, variety and precision of the data used to constrain 
the model. 

The clearest physical indication from the results pre- 
sented in this paper is that our current treatment of star 
formation, although similar to that used both in other phe- 
nomenological models and in direct simulations of galaxy 
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formation, is significantly in error, producing overly efficient 
star formation at early times and in small galaxies. We have 
tried simple modifications of these recipes but have not so far 
identified one which leads to substantially improved results. 
A better astrophysical understanding of large-scale star for- 
mation is probably required. 
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APPENDIX A 

As discussed in Sec. 3.6, we have modified our previous treat- 
ment of the transition between central and satellite status 
when galaxies fall into larger systems. As long as the subhalo 
associated with a galaxy remains outside the virial radius of 
its FOF group, we now continue to treat that galaxy as an 
independent central object. Thus galaxies effectively become 
satellites only when they fall within R v ir- This reduces the 
number of satellites from the point of view of our galaxy for- 
mation modelling and it increases the fraction of satellites 
which are type 2 or "orphan" systems with no associated 
subhalo. (This is because the orphans almost all lie within 
-Rvir.) Here we illustrate the change in the effective number 
of satellites in our two simulations as a function of stellar 
mass. 

Fig. 24 shows the fraction of all galaxies at each stel- 
lar mass which are satellite systems of various types. Red 
curves refer to the MS-II and are plotted down to a stel- 
lar mass of 1O 7 M0, while blue curves refer to the MS and 
stop at its resolution limit, M* ~ 1O 9 ' 5 M0. For each sim- 
ulation the solid curve gives the fraction of galaxies which 
are centred on non-dominant subhalos of their FOF groups, 
the dashed curve gives the fraction which are in addition 
within -Rvir, and the upper and lower dotted curves give the 
fractions which are orphans within FOF groups and within 
-Rvir, respectively. In our previous work (e.g. DLB07) the 
galaxies corresponding to the solid and upper dotted curves 
were treated as satellites when modelling their evolution. In 
the current paper it is rather the galaxies corresponding to 
the dashed and lower dotted curves which are treated as 
satellites; the galaxies corresponding to the difference be- 
tween the solid and dashed curves continue to be treated 
as centrals. Thus the effective satellite fraction is smaller in 
this paper than in our previous work. Notice that while the 
improved resolution of the MS-II does decrease the number 
of orphan galaxies in comparison to the MS, these remain a 
significant population, even at relatively high stellar mass. 
Notice also that above the mass limit of the MS, the total 
fraction of galaxies which are satellites agrees well between 
the two simulations, demonstrating that our treatment of 
orphans in the MS is indeed appropriate, as also concluded 
earlier when discussing Fig. 14. 
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Figure 24. The fraction of all galaxies which arc satellites of 
various types as a function of stellar mass. Blue curves are results 
for the MS and red curves for the MS-II. For each simulation and 
at each stellar mass, a solid curve gives the fraction centred on a 
non-dominant, satellite subhalo (type 1 galaxies), a dashed curve 
gives the fraction which are in addition within -R v i r of halo centre, 
and dotted curves give the fraction with no remaining associated 
subhalo (type 2 or "orphan" satellites; the upper curve refers 
to orphans within the FOF group while the lower only counts 
orphans within R v ir)- Note that in both simulations a substantial 
fraction of the type l's are actually outside il v ir and so continue 
to be treated as central galaxies by our modified prescriptions. 
Note also that while improved resolution reduces the number of 
orphans in the MS-II, they remain a significant population even 
at relatively large stellar mass. 
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